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ANALYSIS  OF  MODULATED  OSCILLATOR  DATA 


1.  INTRODUCTION 

The  Interior  Ballistics  Division  of  the  Ballistic 
Research  Laboratory  has  conducted  a  firing  experiment  with  a 
fully  instrumented  cannon  in  order  to  establish  and  maintain 

a  technology  base  for  the  study  of  projectile-gun  dynamics  . 

Of  particular  importance  is  the  investigation  of  problems  asso¬ 
ciated  with  the  erratic  flight  behavior  of  certain  projectiles. 

The  investigation  of  projectile  in-bore  and  launch 
dynamics  involved  a  number  of  firings  of  an  M-68  105mm  tank 
cannon  during  the  fall  and  winter  of  1977.  The  instrumentation 
of  this  cannon  included  an  inductance  coil  mounted  at  the  muzzle 
face  and  forming  part  of  a  Hartly  radio  frequency  oscillator. 

This  coil  and  oscillator  forms  the  sensor  portion  of  a  measure¬ 
ment  system  which  is  intended  to  provide  an  accurate  muzzle  exit 
time,  monitor  projectile  integrity,  determine  the  projectile 
velocity,  and  trigger  auxiliary  equipment.  Determining  the 
precise  time  of  projectile  exit  and  the  duration  of  the  passage 
of  the  projectile  through  the  coil  are  the  primary  functions  of 
this  instrument. 

The  complete  muzzle  velocity  measurement  system  consists 
of  three  basic  units:  the  sensor,  the  signal  processor,  and  the 
data  processor.  The  sensor  provides  a  pulse  with  amplitude  pro¬ 
portional  to  the  geometric  shape  of  the  projectile.  The  signal 
processor  amplifies  the  detected  pulse,  filters  and  samples  it, 
and  prepares  it  for  data  processing.  The  data  processor  continues 
with  the  extraction  of  time  and  speed  information  from  the  raw 
data . 


The  sensor  is  a  single  loop  inductance  coil  which  has 
been  etched  on  a  single  clad  printed  circuit  board  and  is  rigidly 
attached  to  the  gun  tube  muzzle.  The  configuration  shown  in 
Figure  1.1,  allows  the  sensor  to  be  mounted  without  drilling  holes 
into  the  gun  tube.  The  sensor  coil  is  part  of  the  tuned  circuit 
of  the  Hartly  oscillator  and  therefore  directly  affects  both  the 
amplitude  and  frequency  of  oscillation. 

An  expanded  view  of  the  sensor  coil  arrangement  is 
shown  in  Figure  1.2.  Its  components  are  the  guard  ring,  sensor 
coil,  and  teflon  seal.  The  inductance  coil  is  insulated  and 
separated  from  the  muzzle  face  by  the  thickness  of  the  board 


R.  K.  Coder,  J.  0.  Pilcher,  "Nondestructive  Test  Method  to  Establish  the 
Performance  of  Projectile-Gun  Systems Proceedings  of  the  26-th  Defense 
Conference  on  Nondestructive  Testing,  15  to  17  November  1977 ,  Seattle , 
Washington;  Published  by  U.S.  Army  Materials  and  Mechanics  Research 
Center ,  Watertown,  MA  (pages  272-503). 
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OSCILLATOR 


Figure  1.1  Cross-Sectional  View  of  Sensor  Coil 
Assembly  Rigidly  Attached  to  M-68 
Gun  Tube  Muzzle 


GUARD  RING  TEFLON  SEAL 


Figure  1.2  Sensor  Coil  Expanded  View 


from  which  it  is  made.  A  second  ring  consisting  of  an  unetched 
single  clad  circuit  board  forms  a  shield  (guard  ring)  and  is 
mounted  in  front  of  the  sensor  coil.  Both  pieces  of  circuit 
board  are  affixed  with  resin  to  the  front  face  of  the  collar 
and  securely  held  to  the  gun  tube  muzzle  by  the  mounting  collar. 
To  prevent  high  pressure  gas  from  leaking  between  the  muzzle  face 
and  the  sensor  coil  ring,  a  teflon  seal  had  to  be  inserted.  The 
rest  of  the  oscillator  circuit  is  housed  in  a  small  box  mounted 
on  the  collar  assembly. 

Electronically,  the  oscillator  is  a  conventional  Hartly 
circuit.  The  utilization  of  this  oscillator  for  monitoring  the 
projectile  exit  from  the  gun  muzzle  is  based  on  the  recording  and 
evaluation  of  changes  in  the  electric  current  of  the  oscillator 
which  are  caused  by  the  interaction  of  the  electromagnetic  field 
of  the  inductance  coil  with  the  transient  in— bore  environment 
consisting  of  the  moving  projectile  and  the  gas  flow  which 
surrounds  it. 


Ideally,  the  oscillator  coil  would  have  a  constant 
inductance  and  distributed  capacitance  prior  to  the  arrival  of 
the  projectile,  and  then  would  show  sharp  changes  due  to  the 
magnetic  and  dielectric  properties  of  the  projectile  material. 
Thus  the  desired  time  information  could  be  easily  and  accurately 
extracted  by  simply  comparing  the  sampled  values  on  the  tape 
with  a  fixed  threshold.  More  detailed  analysis  would  yield  a 
profile  of  the  projectile  cross  section  as  a  time  function,  or 
the  relative  positions  of  salient  features  of  the  projectile. 

Since  the  amplitude  of  the  oscillator  current  is 
modulated  as  a  function  of  the  contour  of  the  metallic  shell 
of  the  projectile,  any  discontinuity  in  the  projectile  contour 
can  be  chosen  as  a  geometric  reference  point.  By  correlating 
the  passage  time  of  the  reference  points  with  their  spatial 
displacements,  one  can  determine  the  muzzle  exit  velocitv  of 
the  projectile. 


practice,  the  coil  and  oscillator  sense  the  approach 
of  the  projectile  before  it  arrives  at  the  muzzle,  thus  "softeninq" 
the  change  in  oscillator  output  at  the  time  of  exit.  The  induc¬ 
tance  and  capacitance  of  the  coil  are  affected  by  a  number  of 

phenomena  other  than  the  passage  of  the  projectile  2 ^  .  Gases 
within  the  coil  prior  to  and  following  the  projectile  exit  have 
an  effect  on  the  coil  electrical  properties.  Furthermore,  the 
dimensions  of  the  coil  and  its  proximity  to  the  muzzle  structure 


R.  K.  Loder,  J.  Q.  Schmidt,  "Radio  Frequency  Oscillator  Technique  for 
Monitoring  Velocity  and  Structural  Integrity  of  Projectiles  from  Their 
Ex-it  from  the  Muzzle  »  ARBRL-MR-0 31 00,  USA  ARRADCOM,  Ballistic  Research 
Laboratory ,  Aberdeen  Proving  Ground,  MD,  April  1981.  ( AD  A100725) 
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vary  slightly  during  firing  due  to  vibration  of  the  tube  and  coil 
mounting,  and  thus  the  vibration  is  translated  into  a  time-varying 
inductance  change.  As  it  is  always  the  case  with  electronic 
devices,  the  circuitry  associated  with  the  coil  and  oscillator 
adds  a  random  fluctuation,  or  instrumentation  noise,  component 
to  the  signal.  Generally,  the  vibration  produces  a  quasiperiodic 
component  in  the  recorded  signal,  while  the  instrumentation  noise 
is  nearly  random.  In  actuality,  the  vibration  period  is  not 
constant,  and  the  noise  does  not  necessarily  have  a  Gaussian 
amplitude  distribution  or  flat  power  spectrum.  Thus  very 
simple  filtering  and  smoothing  are  not  sufficient  to  remove 
these  corrupting  components. 

The  instrumentation  noise,  gas,  and  vibration  effects 
on  the  signal  tend  to  even  further  obscure  the  points  where  the 
projectile  arrives  at  and  departs  from  the  coil.  As  an  additional 
problem,  when  severe  saturation  of  the  electronics  occurs  as  a 
result  of  the  signal  exceeding  the  dynamic  range  of  the  circuitry, 
the  relatively  slow  decaying  transients  produced  distort  the 
subsequent  waveform  to  the  extent  that  the  data  are  useless.  In 
some  cases,  exceeding  the  dynamic  range  for  very  brief  periods 
causes  only  minor  clipping  of  the  waveform,  but  even  in  these 
cases  some  data  are  lost. 

While  the  gas  and  vibration  effects  are  a  problem  when 
determination  of  exit  time  and  duration  are  to  be  made,  these 
same  effects  permit  the  recovery  of  some  secondary  data  from  the 
oscillator  waveforms.  The  secondary  data  include  measurement  of 
radial  and/or  longitudinal  gun  tube  vibration,  and  the  location 
of  the  shock  front. 

The  most  important  of  the  coil  loading  phenomena,  for 
the  purpose  here,  is  the  electromagnetic  field  induced  in  the 
projectile  surface  by  the  electromagnetic  field  of  the  oscillator. 
The  physical  process  which  underlies  this  major  loading  source  is 
described  as  follows:  The  varying  electromagnetic  field  radiated 
from  the  active  inductance  coil  produces  localized  eddy  currents 
in  the  metallic  surface  of  the  projectile  when  it  moves  through 
the  coil.  The  eddy  currents,  in  turn,  generate  an  electromagnetic 
field  which  couples  back  into  the  active  inductance  coil.  The 
closer  the  metallic  shell  of  the  projectile  is  to  the  inductance 
coil,  the  stronger  is  the  induced  secondary  field  and  the  effect 
of  its  recoupling  to  the  inductance  coil.  The  variation  in  the 
oscillator  impedance  causes  an  increase  or  decrease  of  the 
operating  frequency  of  the  oscillator  which  can  be  recorded. 

This  frequency  modulation  was  the  basis  for  a  measurement  device 

f  31 

developed  by  Mr.  G.  Schultze  at  ISL  .  However,  with  the  proper 


r  3i 

G.  Schultze,  "Telemetrische  V  - Messungen ";  Colloquium  on  Ballistic 
Measurement  Techniques;  held  at  the  Institut  Franco -Allemand  De 
Recherches  De  Saint-Louis,  6  and  7  June  1973;  ISL  Report  14/73. 
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selection  of  oscillator  circuit  parameters,  the  change  in  im¬ 
pedance  can  be  kept  such  that  the  frequency  shift  is  minimized 
snd  the  amplitude  modulation  is  maximized  for  obtaining  the 
signature  of  the  projectile  exit.  This  technique  of  amplitude 

modulation  was  selected  by  Mr.  J.  Q.  Schmidt for  the  M-68 
instrumentation,  and  is  the  technique  for  which  data  are 
evaluated  by  the  program  described  in  this  report. 


Falcon  Research  and  Development  Company  has  generated 
a  computer  code,  designated  as"MOdulated  Oscillator  Evaluation", 
(MOOSE),  which  uses  the  2.5  microsecond  interval  samples  from 
the  detected  oscillator  signal  as  its  input,  and  performs  the 
following  operations: 


2 

3, 


4, 

5 

6, 


Classification  of  the  recorded  time  series  according 
to  the  quantity  of  data  lost  due  to  saturation: 


a. 

b. 

c. 


No  data  loss. 

Slight  loss  due  to  clipping  during  pass-through, 
Major  section  of  the  history  lost  following 
saturation  during  pass-through. 


Determination  of  the  beginning  and  end  of  pass-through. 
Division  of  the  time  series  into  three  regimes: 

a.  Prior  to  projectile  exit;  Regime  I. 

b.  During  projectile  pass-through;  Regime  II. 

c.  Following  pass-through;  Regime  III. 

Data  smoothing  by  local  mean  or  local  median. 

Determination  of  the  power  spectra  for  Regimes  I  and  III 

Description  of  the  signal  in  Regimes  I  and  III  as  a 
sum  of  three  major  components: 


a , 

b, 


c , 

d, 


An  aperiodic  function  which  is  modeled  as  a  4th 
degree  polynomial. 

A  periodic  component  which  is  modeled  as  the  sum 
of  ten  sinusoidal  functions,  harmonically  related 
but  not  necessarily  by  consecutive  harmonics. 
Random  noise. 

Identification  of  shock— like  waveforms. 


[4] 

J.  Q.  Schmidt ,  "A  Radio  Frequency  Oscillator  Technique  for  Measuring 
Prog ectile  Muzzle  Velocity USARRADCOM ,  Ballistic  Research  Laboratory , 
Aberdeen  Proving  Ground,  MD ;  Technical  Report  ARBRL-TR-02158  April 
1979.  (AD  B0S8926L) 
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Note:  As  an  option,  the  first  regime  may  be  divided  into  two 
parts,  with  the  division  point  being  that  point  where  the 
periodic  component  changes  its  fundamental  frequency.  The 
two  parts  of  the  first  regime  are  then  analyzed  in  the 
manner  described  for  the  complete  regime. 

7.  Determination  of  the  polynomial  coefficients  by  an 
orthogonalized  polynomial  technique. 

8.  Determination  of  the  ten  sinusoidal  frequencies  by 
choosing  the  frequencies  for  which  the  power  spectrum 
shows  the  greatest  energy. 

9.  Checking  the  random  component  for  nonperiodicity  and 
"whiteness"  by  power  spectrum  analysis,  and  for  the 
presence  of  a  Gaussian  amplitude  distribution  by  a 
Chi  Square  (y2)  test. 

10.  Determination,  for  Regime  II,  of  the  following: 

a.  Extrema. 

b.  Power  Spectra. 

The  program  was  prepared  in  FORTRAN  and  is  intended 
for  use  at  the  BRL  computer  facility  at  Aberdeen  Proving  Ground. 
In  the  preparation  of  the  program,  extensive  use  was  made  of 
the  IMSL  subroutines  in  order  to  minimize  coding  time.  The 
mathematical  analysis  performed  by  MOOSE,  and  the  computer  code 
itself,  are  fully  described  in  the  following  sections. 
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2. 


ELEMENTS  OF  ANALYSIS 


2.1  Data  Description. 

Measurements  were  taken  by  the  induction  coil  and 
oscillator  for  a  total  of  five  milliseconds,  beginning  approx¬ 
imately  at  the  time  of  ignition.  The  time  interval  between  data 
points  is  2  5  microseconds.  The  data  are  in  units  proportional 
to  the  demodulated  signal  voltage  and  have  values  ranging  between 
approximately  1800  and  -2300  units,  although  higher  values  did 
occur.  When  the  signal  exceeded  these  limits  for  more  than  25 

microseconds,  saturation  of  the  instrumentation  led  to  a  loss 
of  subsequent  data. 


During  the  time  period  in  which  the  induction  coil 
operates,  three  stages  of  data  are  represented.  The  first  is 
just  after  ignition,  while  the  projectile  is  still  completely 
iLtJev,tUbe‘,  during  the  second  stage,  the  projectile  is  exiting 
the  tube,  and  the  third  is  after  the  projectile  has  completely 
emerged  from  the  tube.  The  signal  shape  for  each  of  these  stages 
is  distinct,  although,  for  the  various  reasons  discussed  in  the 
introduction  the  divisions  between  them  are  not  always  well 
defined.  For  the  purposes  of  analysis,  each  of  the  physical 
stages  of  the  firing  corresponds  to  a  regime  in  the  data  set. 

The  first  regime  begins  with  a  noisy,  relatively  cons¬ 
tant  signal  at  or  below  the  zero  level.  In  some  cases  a  strong 
periodic  signal  is  introduced  during  regime  one,  generally 
accompanied  by  a  shock  signal.  Y 

The  second  regime,  using  the  proof  slug  data  for 
illustration  purposes,  is  marked  by  a  sharp  rise  in  the  level  of 
e  signal,  well  m  excess  of  the  average  noise  level,  as  the 
projectile  enters  the  induction  coil.  The  signal  peaks  and 
ecreases  rapidly  by  approximately  one  third.  This  spike  is 
he  most  prommant  feature  in  the  data  set,  although  it  does  not 
necessarily  represent  the  maximum  signal  value.  Following  the 
signal  decrease,  there  is  a  period  of  constant  value  followed  bv 

he“cSeTller  ^lke-.,The  si9nal  immediately  plummets  as 
.projectile  leaves  the  coil,  perhaps  evincing  one  or  two 

spikes  during  the  descent,  depending  on  the  cross  section  of  the 
Lie  The  signal,  after  reaching  its  minimum? moves 
2^ e^tly  into  the  third  regime.  This  minimum  is  generally  not 
far  from  the  baseline  signal  level  in  regime  one.  Conventional 
types  of  projectiles  show  a  different  form  for  the  second  regime. 
The  noisy  interval  between  the  two  peaks  is  reduced  to  a  point 
or  two  and  its  level  is  very  near  that  of  the  two  pea?s?  ??n 
addition,  no  spikes  are  exhibited  as  the  signal  attains  its 
minimum  at  the  end  of  the  regime. 

rr  Finally,  after  the  projectile  has  completely  emerqed 

from  the  gun  barrel,  the  induction  coil  measures  the  signtll 
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still  being  produced  by  gases  and  vibrations.  This  is  in 
Regime  III  of  the  data  set.  There  are  three  basic  components 
to  this  interval.  The  first  is  the  baseline  which  rises  from 
the  minimum  at  the  last  point  in  the  second  regime  and  even¬ 
tually  either  decreases  or  levels  off.  Superimposed  on  this 
baseline  is  a  well-defined  periodic  signal.  Instrumentation 
noise  occurs  throughout  Regime  III,  possibly  punctuated  by 
shocks. 


Representative  data  sets  are  plotted  in  Figure  2.1, 
showing  typical  characteristics  of  Regimes  I,  II,  and  III,  and 
the  differences  in  the  signals  of  two  types  of  projectiles. 

2.2  Requirements  of  Analysis. 

As  discussed  in  the  previous  section,  there  are 
several  distinct  signal  components  present  in  the  induction 
coil  data,  each  representing  some  physical  process  in  the  gun 
firing.  It  is  assumed  that  these  components  add  linearly, 
although  in  some  cases,  such  as  vibration  of  the  coil  windings, 
some  multiplicative  combinations  may  occur.  The  requirements 
for  separating  the  signal  components  and  the  implications  of 
these  requirements  are  discussed  in  the  following  paragraphs. 

The  first  requirement  is  that  the  computer  program 
automatically  identify  any  data  sets  in  which  the  instrumentation 
became  overloaded  or  saturated.  A  decision  is  then  necessary 
as  to  the  extent  of  the  analysis  which  is  possible  for  that  data 
set,  and  the  amount  of  data  which  is  unaffected  by  the  saturation. 

After  the  identification  of  incomplete  or  "damaged" 
data  sets,  the  next  function  of  the  analysis  is  to  determine 
the  duration  of  Regime  II,  which  is  the  time  taken  by  the  pro¬ 
jectile  to  pass  the  induction  coil.  Identification  of  the 
epochs  defining  the  end  points  of  Regime  II  also,  of  course, 
determines  the  end  of  Regime  I  and  beginning  of  Regime  III. 

Because  knowledge  of  the  epoch  locations  is  critical  to  the 
calculation  of  muzzle  velocity,  an  estimate  should  be  made  of 
the  error  expected  in  determining  their  locations. 

The  next  priorities  are  the  analyses  of  the  first  and 
third  regime.  Regime  I  may  be  divided  into  two  subregimes  char¬ 
acterized  by  a  pronounced  shift  in  the  fundamental  frequency  of 
the  periodic  component,  but  the  same  general  process  is  used  on 
each  data  interval.  This  analysis  separates  the  four  components 
of  the  signal.  The  first  component  is  the  trend,  or  baseline 
of  the  interval.  The  second  is  the  periodic  signal.  Last  are 
the  noise  and  the  occasional  shocks.  Equations  are  to  be  found 
to  describe  the  baseline  and  the  periodic  component.  The  noise 
is  described  in  terms  of  its  spectrum  and  distribution  of  the 
amplitudes,  rather  than  by  an  equation.  Points  that  may  repre¬ 
sent  shocks  are  indicated. 
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Nonrandom  trends  or  periodicities  in  the  "noise"  may 
reveal  inadequacies  in  the  derivation  of  the  baseline  and  peri¬ 
odic  function  of  a  regime. 

Finally,  there  may  be  two  different  portions  to 
Regime  I.  The  transit  of  the  shock  front  of  the  compressed 
air  in  front  of  the  projectile  may  manifest  itself  as  a  shock 
in  the  signal.  When  this  occurs,  the  signal  waveform  succeeding 
the  shock  shows  more  baseline  variation  and  periodicities  than 
the  anteceding  one .  To  adequately  model  Regime  I ,  it  may  be 
advantageous  to  apportion  it  into  two  parts,  with  each  analyzed 
separately.  The  requirement  for  accuracy  in  detecting  this 
epoch  is  less  critical  since  each  subregime  can  be  adequately 
modeled  and  analyzed  even  when  the  chosen  division  point  is 
slightly  in  error. 

Table  2-1  summarizes  all  of  the  requirements  of  the 
analysis  for  the  data. 


TABLE  2-1  REQUIREMENTS  FOR  THE  ANALYSIS 
OF  THE  OSCILLATOR  SIGNALS. 


1.  Classification  into  good  or  "damaged" 
data  sets. 

2.  Length  of  time  in  Regime  II. 

3.  Periodicities  of  Regimes  I  and  III. 

4.  Description  of  noise  of  Regimes  I 
and  III. 

5.  Baseline  of  signal  of  Regimes  I 
and  III. 

6.  Relative  location  of  shocks. 

7.  Accuracy  of  epoch  identification. 

8.  Location  of  preliminary  shock  in 
Regime  I . 

9 .  Power  spectrum  and  extrema  of 
Regime  II. 
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2 . 3  APPROACH . 


In  this  section,  the  major  assumptions  underlying  the 
approach  to  the  four  major  phases  of  the  analysis  of  the  oscill¬ 
ator  signals  will  be  discussed.  These  four  phases  include 
identification  of  incomplete  or  damaged  data  records,  location 
of  the  transition  points  between  the  three  data  regimes,  analysis 
of  the  data  in  each  regime,  and  analysis  of  the  residual  noise. 

Identification  of  damaged  data  records  is  made  by  the 
algorithm  shown  in  Figure  2.2.  This  algorithm  essentially  looks 
for  indications  that  the  recording  electronics  became  saturated 
during  Regime  II,  that  is,  during  the  passage  of  the  projectile 
through  the  induction  coil.  It  is  assumed  that  the  saturation 
signal  will  have  the  characteristics  which  this  routine  detects, 
in  particular  that  successive  data  points  will  have  the  same 
amplitude  during  the  saturation  portion  of  the  record.  A 
further  assumption  made  for  this  case  is  that  short-term  satu¬ 
ration,  less  than  10  points  or  25  microseconds,  or  "clipping" 
will  have  relatively  little  effect  on  the  signal  following  the 
return  of  the  recorded  values  into  the  linear  range  of  the 
instrumentation.  Furthermore,  it  is  assumed  that  heavy  satura¬ 
tion,  as  detected  by  this  algorithm,  will  cause  transients  in 
the  instrumentation  output  which  render  the  remaining  portion 
of  the  record  useless,  and  the  MOOSE  program  automatically 
analyzes  only  the  first  regime. 

The  second  phase  of  the  analysis  is  the  identification 
of  the  transition  points  between  Regimes  I  and  II  and  between 
Regimes  II  and  III.  These  transition  points,  or  "epochs",  are 
produced  when  the  projectile  enters  and  leaves  the  induction 
coil,  and  thus  represent  transitions  between  three  physically 
distinct  portions  of  the  data.  Detection  of  the  epochs  is  the 
essential  operation  in  determining  the  projectile  exit  velocity, 
and  for  providing  key  timing  information  for  analyzing  the  data 
from  other  instrumentation.  It  must  be  assumed  that  the  epochs 
can  be  determined  with  the  accuracy  needed  for  these  purposes 
by  either  operator-assisted  interactive  means  or  by  fully 
automatic  techniques,  and  furthermore  it  is  assumed  that  no 
data  other  than  that  provided  on  the  oscillator  coil  record  is 
needed  for  the  determination. 

Research  in  signal  analysis  has  been  concerned  with 
epoch  detection  of  this  type  for  many  years.  In  general,  the 
approach  has  been  to  fit  the  data  with  a  function,  or  sum  of 
functions,  and  then  for  each  point  in  the  record  extrapolate 
the  data  both  forward  and  backward.  An  epoch  caused  when  a 
new  signal  is  added  to  the  existing  one  will  be  indicated  by 
a  large  discrepancy  between  the  values  obtained  by  forward 
and  backward  extrapolation.  This  approach  works  well  when 
the  signals  are  sums  of  exponentials  and  the  noise  is  Gaussian. 
This  technique  was  attempted  with  the  gun  data  but  was  found 
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Data  Record 


Figure  2.2  Data  Classification  Algorithm 


to  give  erratic  results.  A  technique  based  on  measurement  of 
the  curvature  of  the  signal  in  the  region  where  the  epoch  was 
expected  gave  even  more  disappointing  results.  The  approach 
finally  chosen  for  this  determination  is  discussed  in  detail 
in  the  next  section  of  this  report. 

The  third  phase  of  the  analysis  is  the  actual  extrac¬ 
tion  of  the  additive  components  of  the  signals  in  Regimes  I  and 
III.  The  signals  are  broken  into  polynomial,  periodic,  and 
random  components  by  straightforward  curve  fitting  and  Fourier 
techniques.  The  major  assumption  in  this  phase  is  that  the 
polynomial  component,  which  is  a  measure  of  the  shift  of  the 
data  baseline  with  time,  and  the  periodic  function,  which  is  a 
measure  of  the  vibration  of  the  tube  during  firing,  are  additive 
in  their  effects  on  the  signal  from  the  induction  coil.  This 
implies  linearity  in  measurement  and  recording  equipment  as  well 
as  independence  of  the  physical  effects.  The  actual  approach 
used  in  subtracting  out  the  components  of  the  signal  is  described 
in  the  next  section  and  represents  the  major  phase  of  the  analy¬ 
sis  presented  here. 

The  final  phase  of  the  analysis  is  the  determination 
of  the  characteristics  of  the  random  "noise"  component  of  the 
data.  The  residual  signal  after  both  the  polynomial  and  periodic 
components  have  been  successfully  subtracted  out,  is  assumed 
to  be  purely  random,  although  its  amplitudes  are  not  assumed 
necessarily  to  be  Gaussian,  nor  are  the  amplitudes  necessarily 
uncorrelated  in  the  short  term  sense.  Stated  another  way,  it 
is  not  assumed  that  the  power  spectrum  of  the  residual  noise 
will  be  flat  or  "white". 


19 


3. 


METHODOLOGY 


3 . 1  PROCEDURE ♦ 

The  first  step  in  the  analysis  is  the  identification 
of  the  two  epochs  defining  the  beginning  and  end  of  Regime  II. 

As  discussed  in  Section  2.1,  the  second  regime  has  certain  well- 
defined  features,  and  these  can  be  exploited  in  determining  the 
beginning  and  ending  points.  This  regime  begins  with  a  sharp 
rise  in  the  signal  level  to  a  value  greatly  in  excess  of  the 
noise  level.  The  signal  peaks,  then  falls  slightly  before 
rising  to  a  second  somewhat  lower  peak.  Following  this  second 
peak,  the  signal  decreases  rapidly  and  may  contain  one  or  two 
short  duration  spikes.  The  signal  reaches  a  minimum  at  the  end 
of  Regime  II,  and  moves  into  Regime  III  which  is  essentially  a 
noisy  periodic  function  on  a  baseline  which  initially  is  in¬ 
creasing.  The  M392  projectiles  produce  a  similar  pattern, 
although  the  noise  between  the  two  peaks  is  less  and  there  is 
an  absence  of  spikes  during  the  signal  decrease. 

A  decision  must  be  made  as  to  whether  the  epochs  should 
be  determined  solely  from  this  waveform,  making  no  use  of  any 
knowledge  of  the  underlying  physical  phenomena,  or  whether  an 
attempt  should  be  made  to  deduce  the  expected  signal  from  the 
physical  apparatus  and  use  this  expected  signal  as  a  basis 
for  determining  Regime  II.  Since  a  better  understanding  of  the 
physical  phenomena  is  the  purpose  for  which  the  experiment  was 
performed,  it  was  decided  that  predicting  the  waveform  in  order 
to  determine  the  regimes  was  an  unjustified  bootstrapping,  and 
the  determination  of  epochs  was  therefore  done  solely  on  the 
basis  of  the  observed  waveform.  The  assumption  was  made  that 
the  waveform  consists  of  two  major  additive  components:  a 
signal  consisting  of  noise,  periodic  functions,  and  a  growing 
baseline,  to  which  is  added  the  Regime  II  signal.  There  are 
two  characteristics  in  this  overall  pattern  which  are  useful  in 
identifying  the  epochs.  The  first  is  the  high  signal-to-noise 
ratio  on  the  increasing  and  decreasing  waveforms  near  the 
boundaries  of  Regime  II.  The  second  is  the  high  slope  of  the 
signal  near  its  first  peak.  Obscuring  the  epochs,  however,  are 
the  noise,  which  is  especially  evident  in  Regimes  I  and  III,  the 
shock  signals  at  the  end  of  Regime  II,  and  the  similarities 
between  the  periodic  waveform  at  the  end  of  Regime  I  and  the 
(gradually)  increasing  slope  in  the  signal  at  the  beginning  of 
Regime  II.  Any  procedure  which  is  used  to  detect  the  epochs 
must  take  these  into  account  while  using  the  identifying 
features  of  the  Regime  II  epochs. 

Several  approaches  to  epoch  determination  were 
attempted  and  rejected  before  the  procedure  described  here  was 
decided  upon.  The  first  and  simplest  was  to  measure  the 
curvature  of  the  signal  near  the  expected  epoch  points,  and 
choose  that  point  where  the  curvature  was  greatest.  With 
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actual  signal  data,  this  method  was  found  to  be  unstable  and 
generally  converged  to  a  point  which  was  intuitively  unsatisfying. 
A  second  approach  was  to  fit  polynomials  to  the  data  well  within 
Regime  II,  and  well  into  Regime  I,  and  then  extrapolate  forward 
from  Regime  I  and  simultaneously  extrapolate  backwards  from 
Regime  II,  in  the  hope  that  the  largest  discrepancy  between 
the  extrapolated  values  would  occur  at  the  true  epoch.  Results 
with  this  approach  were  inconsistent,  and  the  computer  time 
required  was  considered  excessive.  Methods  based  directly  on 
differentiating  the  signal  were  unsuitable  because  the  differen¬ 
tiating  had  a  markedly  deleterious  effect  on  the  signal-to-noise 
ratio,  and  thus  further  obscured  the  epoch  which  was  to  be 
detected . 


The  algorithm  which  was  actually  used  to  detect  the 
epochs  is  shown  in  Figure  3.1.  As  a  reference  point  from  which 
the  first  and  second  epochs  are  located,  the  point  at  which  the 
signal  has  its  maximum  first  derivative  (maximum  slope)  was  used, 
giving  a  point  just  before  the  peak.  (This  peak  is  not  necess¬ 
arily  the  maximum  signal  achieved  during  the  shot.)  The  first 
epoch  is  found  by  a  two-step  process.  Beginning  with  the 
reference  point,  and  considering  only  the  preceding  two  hundred 
points,  points  are  successively  examined  until  one  is  found  which 
is  less  than  one  of  the  points  preceding  it  in  that  interval. 

This  is  the  initial  estimate  of  the  first  epoch.  The  final 
estimate  of  the  first  epoch  is  determined  to  be  the  first  point 
in  the  interval;  bounded  by  the  initial  estimate,  that  begins  a 
monotonically  increasing  sequence  up  to  the  estimate.  This 
epoch  is  accurate  to  within  about  three  points  in  either  direction. 

The  second  epoch  is  found  by  finding  the  first  occurrence 
of  the  minimum  whose  preceding  three  points  are  monotonically 
decreasing,  within  a  two  hundred  point  interval  beginning  at  the 
point  of  the  maximum  first  derivative.  Accuracy  is  estimated  to 
be  within  one  or  two  points. 

Before  analysis  can  begin  on  the  first  regime,  a 
secondary  epoch  detection  may  be  necessary  within  this  subset 
of  the  data.  This  is  required  when  there  is  a  significant  change 
in  the  fundamental  frequency  of  the  periodic  component  of  the 
signal.  The  user  decides  whether  or  not  such  a  change  occurs 
and,  if  so,  the  point  of  division  by  subjective  analysis  of  the 
data  set  and/or  the  examination  of  the  power  spectra  of  various 
subregimes.  If  this  procedure  is  necessary,  there  will  be  three 
subsets  of  data,  instead  of  two  which  must  be  analyzed  separately. 

The  analysis  for  the  third  regime  and  for  the  first 
regime  or  its  two  sections  are  identical.  The  object  is  to 
separate  the  periodic  function,  its  baseline,  and  the  shocks  and 
noise  along  the  curve.  The  baseline  is  assumed  to  be  a  low 
degree  polynomial .  A  fourth  degree  polynomial  was  found  to  be 
the  lowest  degree  to  give  a  reasonable  fit  to  the  curve.  The 
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Preliminary  Data  Record 


Figure  3.1  Epoch  Detection  Algorithm 


,  NPT 
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technique  used  for  fitting  used  a  least  squares  criterion  with 
orthogonalized  polynomials.  The  resulting  polynomial  was  sub¬ 
tracted  from  the  original  data,  leaving  the  periodic  function 
and  noise.  A  power  spectrum  generated  for  this  remainder  curve 
identified  the  ten  frequencies  with  the  most  power.  These  were 
then  used  to  create  a  periodic  function.  The  value  of  ten  is 
somewhat  arbitrary  but  is  believed  to  give  a  good  balance  be¬ 
tween  accuracy  and  simplicity.  The  actual  residual  error 
cannot  be  estimated  by  elementary  Fourier  techniques  because 
of  the  presence  of  the  noise.  This  process  is  described  in 
more  detail  in  Section  3.2,  Basic  Analysis  Techniques.  When 
the  periodic  functions  and  the  polynomial  have  been  subtracted 
from  the  original  curve,  the  remaining  signal  consists  primarily 
of  noise,  with  some  shocks.  When  the  periodic  function  is  added 
to  the  polynomial ,  the  result  is  a  good  approximation  to  the 
original  noiseless  signal. 

In  the  final  step  of  the  analysis,  the  noise  in  the 
first  and  third  regimes  is  examined.  The  power  spectrum  for  each 
is  generated.  The  noise  is  then  characterized  as  white  noise  or, 
if  it  is  non-white,  its  dominant  frequencies  are  specified.  A 
X  test  on  the  noise  gives  the  mean  and  standard  deviation  of 
the  noise  and  indicates  the  probability  of  the  noise  being 
Gaussian. 


Identification  of  the  shocks  remaining  in  the  noise 
is  done  by  determining  the  points  greater  than  2.5  standard 
deviations  from  the  mean.  These  are  printed  and  points  which 
are  more  than  3  standard  deviations  from  the  mean  are  specially 
noted.  In  addition,  the  number  of  standard  deviations  of  the 
first  derivative  from  its  mean  are  printed  for  those  points. 

It  is  expected  that  a  shock  occurs  when  both  measures  give  high 
deviations  and/or  when  adjacent  points  are  far  from  the  mean  in 
opposite  directions. 

This  concludes  the  analysis  of  the  gun  signal. 

3.2  BASIC  ANALYSIS  TECHNIQUES. 

Most  of  the  analysis  techniques  required  for  the 
inductance  coil  data  are  standard  mathematical  and  statistical 
procedures.  A  few  techniques,  including  those  for  the  detection 
of  damaged  data  and  for  epoch  detection,  are  peculiar  to  this 
particular  problem,  and  these  required  the  development  of 
specialized  computer  subroutines.  In  this  section,  techniques 
and  subroutines  acquired  by  Falcon  for  the  coil  data  analysis 
will  be  discussed.  Not  all  of  the  techniques  described  here 
are  currently  used  in  the  data  analysis,  but  they  have  been 
included  because  of  their  past  or  potential  usefulness  in 
related  work. 
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Three  of  these  techniques  are  taken  directly  from 
the  collection  of  subroutines  written  and  issued  by  the 
International  Mathematical  and  Statistical  Libraries,  Inc. 

This  collection,  known  as  the  "IMSL"  library,  is  available  in 
a  number  of  data  processing  centers  throughout  the  country. 

The  BRL  CYBER  computer  center  currently  subscribes  to  the  IMSL 
library.  Due  to  the  large  data  arrays,  it  was  necessary  to 
create  a  special  library  of  IMSL  routines  for  MOOSE,  in  which 
subroutine  storage  arrays  were  increased  and  stored  in  Level  2. 
All  IMSL  routines  needed  by  MOOSE  were  in  this  IMSL 2  library. 

The  first  basic  technique  concerns  time  series  analysis 
in  the  time  and  frequency  domains.  In  particular,  a  power 
spectrum  representation  for  a  given  time  series  is  needed  to 
identify  the  dominant  frequencies  in  the  data  set.  This  is 
accomplished  by  the  IMSL  subroutine  FTFREQ.  For  a  given  time 
series,  this  subroutine  computes  the  mean  and  variance,  the 
autocovariance  function,  and  the  power  spectrum.  The  resolution 
of  this  subroutine  is  250  Hz  when  the  time  series  contains  at 
least  800  points  taken  at  2.5  microsecond  intervals,  decreasing 
to  1  kHz  resolution  when  only  200  points  are  available. 


Given  the  power  spectrum  of  the  periodic  component  of 
the  data,  a  "noise  free"  approximation  of  this  component  was 
synthesized  from  ten  sinusoids.  The  frequencies  for  these 
sinusoids  were  selected  to  be  those  ten  frequencies  for  which 
the  power  spectrum  showed  the  greatest  values.  The  phases  and 
amplitudes  for  the  sinusoidal  components  were  determined  by 
the  following  process: 


x j  =  point  in  original  curve  j  =  1,  NPT  (#  points) 

f  =  desired  frequency  i  =  1,  NFR  (#  frequencies) 

AC.,  BS.  =  coefficients  of  cosine  and  sine  components  of  f. 
i'  l  1 

NPT 

AC.  =  2  l  x.  cos  (2tt  ( j-1)  f  .  At) 

1  j=l  11 

NPT 

BS  .  =  2  l  x.  sin  (2tt  ( j-1)  f  .  At) 

1  j=l  3 

PER.  =  approximation  curve  j  =  1,  NPT 


PER. 

3 


1  NFR 

=  mean  (x)  +  _  £ 


AC^  cos  (2tt  ( j-1)  f ^At)  + 
BS^  sin  (2tt  ( j-1)  f^At) 
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The  technique  for  fitting  a  polynomial  to  the  data 
set,  in  order  to  obtain  a  "noise  free"  approximation  to  the 
baseline,  was  selected  directly  from  the  IMSL  library.  The 
subroutine  chosen  was  RLFOR,  and  uses  orthogonal  polynomials 
in  a  univariate  curvilinear  regression  model  to  fit  a  data 
set.  The  user  specifies  the  maximum  degree  for  the  polynomial, 
the  weight  of  each  point,  and  a  limiting  value  for  the  measure 
of  goodness  of  fit.  This  limiting  value  defines  the  degree 
of  the  polynomial.  The  output  includes  the  coefficients  of 

•  the  fitted  polynomial  and  the  values  of  the  polynomial  function 

for  the  data  set. 

,  The  third  IMSL  subroutine  is  GFIT ,  which  tests  the 

goodness  of  fit  of  a  data  set  to  a  particular  distribution. 

The  x  test  is  performed  using  the  number  of  cells  specified 
by  the  user.  The  distribution  used  for  comparison  is  currently 
the  Gaussian  distribution.  The  pertinent  output  is  the  x2 
statistic  and  the  Q-statistic,  which  is  the  probability  that 
X  is  exceeded  if  the  data  do  not  fit  the  distribution.  In 
other  words,  a  small  Q  implies  that  the  data  fit  the  distri¬ 
bution,  while  a  large  Q  implies  that  they  do  not. 

Before  the  analysis  of  a  data  set  can  begin,  it  must 
be  determined  whether  the  data  are  suitable  for  complete  analysis, 
or  have  been  damaged  through  instrumentation  saturation  to  the 
extent  that  only  partial  analysis  can  be  accomplished.  The 
data  set  classification  checks  for  multiple  occurrences  of 
the  maximum  signal.  There  are  three  classes  of  data  sets: 
no  saturation,  limited  saturation  with  less  than  ten  points 
at  the  maximum  value,  and  heavy  saturation  with  at  least  ten 
points  at  the  maximum  value.  The  subroutine  which  makes  this 
determination  was  programmed  by  Falcon  using  an  algorithm 
supplied  by  BRL . 

Several  data-smoothing  options  exist.  These  include 
the  11-point  moving  average  and  the  11-point  moving  median. 

These  smoothing  techniques  are  not  needed,  nor  are  they  used, 
in  the  current  analysis  since  the  polynomial  and  sinusoidal 
fitting  techniques  inherently  smooth  the  data. 


The  calculations  of  the  first  and  second  derivatives 
are  accomplished  by  means  of  central  differences.  The 
equations  are: 


f ' (*)  =  (8f(x+h)  ~  8f(x-h)  -  f(x+2h)  +  f (x-2h) ) 

f"(x)  =  (16f  (x+h)  +  16f  (x-h)  -  30f  (x)  -  f(x+2h) 

12h 


(1) 

f  (x-2h)  )  (2) 
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In  addition,  routines  exist  which  take  the  mean, 
median,  or  standard  deviation  of  a  set  of  points.  One  routine 
sorts  the  rows  of  two  column  vectors  so  that  the  second  column 
is  in  descending  order. 

Finally,  it  is  possible  to  choose  a  subset  of  the 
data  set  by  specifying  the  beginning  and  ending  times  of  the 
subset  and  the  interval  between  the  points  in  the  data  set. 

3.3  PLOTTING  CAPABILITIES. 

Without  question  the  most  efficient  means  for  pre¬ 
senting  the  results  of  the  analysis  of  a  large  data  set  is 
through  graphs  and  charts.  A  salient  feature  of  the  work 
described  here  is  the  generation  of  plots  at  a  number  of  stages 
of  the  analysis.  For  this  purpose,  the  CALCOMP  Model  780 
digital  plotting  system  was  used.  This  consists  of  a  model 
770  tape  system  and  a  model  763  plotter.  Two  packages  of 
plotting  subroutines  are  used:  the  CALCOMP  plotting  routines 
for  the  BRL ,  CDC  computing  system,  and  the  BRLESC  FORTRAN 

plotting  subroutines  ^ .  With  these  plotting  packages,  graphs 
can  be  plotted  in  a  variety  of  ways.  One  or  more  graphs  can  be 
plotted  on  one  page  using  axes,  borders,  or  grids,  with  either 
linear  or  logarithmic  spacing.  The  points  of  the  graph  can  be 
represented  by  any  one  of  a  large  number  of  symbols  and  can  be 
connected,  if  desired  by  solid  or  broken  lines.  Titles  and 
headings  can  be  included  on  the  plot.  A  major  advantage  of 
this  plotting  system  is  that  the  page  can  be  as  large  as  29 
inches  by  120  feet,  so  that  the  resolution  possible  is  limited 
only  by  the  2.5  microsecond  data  sampling  rate  and  not  by  the 
plotter  itself. 

There  are  currently  four  plotting  options  which  can 
be  utilized.  The  first  is  the  power  spectrum  for  the  first  25 
frequencies,  plotted  on  a  semilog  graph,  with  log10  power 

plotted  against  frequency.  The  second  option  is  a  graph  of  the 
entire  set  of  raw  data,  along  with  the  first  and  second  deri¬ 
vatives  of  the  data.  The  third  option  shows  the  analysis  of  a 
subset  of  the  data,  and  consists  of  three  pages.  The  first 
page  contains  a  plot  of  the  raw  data  in  a  solid  line  and  a 
broken  line  for  the  polynomial  fit.  The  second  page  gives  the 
raw  data  minus  the  polynomial  as  a  solid  line ,  and  the  per¬ 
iodic  fit  as  a  broken  line.  The  third  page  contains  only  the 
noise,  which  is  the  polynomial  and  periodic  fits  subtracted 
from  the  raw  data.  The  fourth  plotting  option  shows  the  entire 
data  set  with  the  noise  subtracted  from  the  first  and  third 
regimes . 


^  Monte  W.  Coleman,  John  V.  Lanahan ,  "BRLESC  FORTRAN  Plotting  Subroutines 
ARDC-TR-6 ,  July  1970,  AD  711893. 
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4.  MOOSE  PROGRAM 

The  MOOSE  (Modulated  Oscillator  Evaluation)  Program 
is  listed  in  Appendix  A.l.  The  program  was  written  in  FORTRAN 
IV,  compatible  with  the  BRL  CYBER  computer  system,  using 
structured  programming  techniques.  The  CALCOMP  plotting 
library  is  required  for  this  program  along  with  IMSL2 ,  a 
special  library  of  IMSL  routines,  which  have  been  converted 
to  level  2  storage.  Appendix  A. 2  lists  the  necessary  IMSL 
routines  which  are  included  in  the  IMSL2  library. 

The  MOOSE  program  uses  level  2  storage  and  organizes 
subroutines  into  different  overlays  for  loading  because  of  the 
size  of  the  program  and  the  storage  required  for  it.  Figure  4.1 
shows  the  overall  flow  of  the  program.  Figure  4.2  shows  the 
organization  of  the  overlays  used  in  the  program.  The  size 
and  storage  requirements  mean  that  the  program  can  only  be 
run  on  mainframe  "Z".  The  program  in  UPDATE  form  as  MOOSEOLDPL 
is  located  on  mainframe  "A" .  The  inductance  coil  data  is  in 
a  file  on  mainframe  "Z"  and  is  attached  as  TAPEl . 

The  job  control  language  for  running  the  MOOSE  pro¬ 
gram  is  shown  in  Appendix  B.l.  Appendix  B.2  gives  a  des¬ 
cription  of  the  input  data  and  Appendix  B.3  shows  an  input 
sample.  The  output  description  of  the  program  is  listed  in 
Appendix  C.l  and  Appendix  C.2  contains  an  output  sample. 
Execution  error  codes  are  listed  in  Appendix  C.3. 
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Figure  4 . 1  MOOSE  Program  Flowchart 


KITTY 
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Figure  4.1  MOOSE  Program  Flowchart  (Continued) 


COMPLT 
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Figure  4.1  MOOSE  Program  Flowchart  (Continued) 


Figure  4.1  MOOSE  Program  Flowchart  (Continued) 
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Figure  4.2  MOOSE  Program  Overlay  Organization 
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APPENDIX  A 

THE  COMPUTER  PROGRAM 
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APPENDIX  A. 1 

TABULATION 
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o  n  n 


c - - - - - 

C  USER  DIAGNOSTICS 

CLOCATION  CODE  MUUE/RE  SUL T  USER  DIAGNOSTIC  MESSAGE 

C- - - -  - 

CSk-SELECT  2  BATCH/FATAL  IMPROPER  SELECTION  MAUE  FOR  PkOCESSING 

CSk-SELECT  2  1A  /INF  OPTION.  CUkkENTLY  ALLOWED  OPTIONS  ARE: 

CSR-SELECT  2  END  -TO  END  PkOCESSING 

CSk-SELECT  2  ANALYSIS  -TO  PERFORM  AN  ANALYSIS 

CSk-SELECT  1  bATCH/FATAL  IMPROPER  SELECTION  MADE  FOR  RUN  MODE. 

CSR-SELECT  1  IA  /FATAL  PROPER  SELECTIONS  ARE: 

CSR-SELECT  1  0  -  FOR  INTERACTIVE  MODE 

CSk-SELECT  1  1  -  FOR  BATCH  MODE 

CSR-SELECT  1  FOR  THE  VARIABLE  IFBACH 

CSk-ANALI  Z  a  BATCH/FATAL  IMPROPER  SELECTION  MAUE  FOk  ANALYSIS 

CSk-ANALIZ  a  IA  /NA  OPTION.  CURRENTLY  ALLOWED  OPTIONS  ARE: 

CSk-ANALI Z  H  0  FOR  EITHER  MODE 

CSk-ANALIZ  9  FOR  THE  VARIABLE  1ANOPT 

CSR-ANALIZ  b  BATCH/NA  IMPROPER  SELECTION  MADE  FOR  ANALYSIS 

CSK-ANALIZ  b  IA  /INF  OPTION.  CURRENTLY  ALLOWED  OPTIONS  ARE: 

CSk-ANALIZ  5  0  -  FOR  EITHER  MODE 

CSR-ANALIZ  5  FOR  THE  VARIABLE  1ANOPT 

OVERLAY  (  A  ,  0 , 0 1 

PROGRAM  MOOSE  ( I NPU T , OUTPUT , T APE b  = 1 NPU T , T APEfe=OUT PUT , T APE1 , TAPE2, 

1  TAPE3,TAPE10,DEBUG=OUTPUT) 

C - - - 

c 

c  *  THIS  PROGRAM  PERFORMS  ONE  OR  MORE  USER  SELECTED  TRANSFORMATIONS 
C  *  ON  GUN-PROJECTILE  DYNAMICS  DATA. 

C 

C  *  FILE  HANDLING  INFORMATION:  TAPES  -  USER  INPUTS  ON  "INPUT" 


C  *  TAPE6  -  RUN  TIME  OUTPUT  SUMMARY 

C  *  T APE  1  -  GUN-PROJECTILE  DATA  FROM  PF 

C  *  TAPE2  -  LOCAL  FILE  TO  HOLD  RUN  INFO 

C  *  TAPE3  -  OUTPUT  FROM  DATA  TRANSFORMATION 

C  *  WHICH  USER  MAY  COPY  TO  OUTPUT 

C  *  OR 

C  *  CATALOG  FOR  RETENTION 

C  *  T APE  1 0  -  USED  FOR  CALCOMP  PLOTTING 

C 

C - FALCON  k&D  JLW  a/19/7 V 


COMMON  /  SSDATA  /  NSSPTS  ,  SOT  I  ME  ,  SUBSTk(SOOO)  ,  SSTIME(5000) 
COMMON  /  SDATA  /  NSTPTS  ,  D  T I M  E  ,  STRING(SOOO)  ,  STIME(bOOO) 
COMMON  /  FREDS  /  FREGEN  (.20)  ,  NFR 

.  COMMON  /REPDAT  /  REPSTR(5000)  ,  I F  REP  ,  NORSTR(50l)0I 
COMMON  /  LABEL  /  KTB(2a)  ,  TST  ART  ,  1 

COMMON  /  DERIVS  /  SlRDOT(bOOO)  ,  STRDD(bOOO)  ,  PLOTOP 
COMMON  /  CLASCM  /  I  CLASS  ,  NBEGBL 

LEVEL  2,  NSSPTSr  SOT  I  ME »  SUBSTk,  SSTIME,  NSTPTS,  DTIME, 

A  STRING,  ST 1ME ,  FREQEN,  NFR, 

B  REPSTR, IFREP,NORSTR,KTB,TSTART,I, 

C  ST  ROOT, ST ROD, PLOTOP, 1  CL  ASS, NBEGBL 

*  WRITE  HEADING  FOR  MOOSE 

WRI TE (b , 900 1 ) 

C 

C  *  START  UP  PROCESSING. 

PROGRAM  LISTING 
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n  n  n  o  o  n  n 


LALL  STAK1 

*  SELECT  UStK  OPTION  FOR  DATA  TRANSFORMATION  DESIRED. 

CALL  ANALIZ 

*  FINISH  RUN 


9000  STOP  "*  NORMAL  TERMINATION  *" 

9001  FORMAT (lHl,///,20X,50(lH*),/,20Xf 1 H* , 52X ,  1 H* , /  ,  20X ,  1H*  , 

*  22X,9hm  0  0  S  E,21X,1H*,/,20X,1H*,52X,1H*,/, 

*  20X,S0(1H*J) 

END 


PROGRAM  LISTING  (Continued) 
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ooo  ooo  o  o  o  o  o  o  o  n  o  o  o 


*  THIS  BLOCK  KATA  SUBPROGRAM  SETS  INITIAL  DATA  VALUES  AT  DEFAULT  LEVEL 

- FALCON  RKD  JLW  a/20/79 

BLOCK  DATA  INITAL 

*  FILE  RELATED  DATA 


*  PROGRAM  CONTROL  VARIABLES 

COMMON  /  CONTRL  /  LU1N  ,  LUOUT  ,  LUDATA  ,  LURUN  ,  LUREPT 

COMMON  /  CUNTRL  /  IFSVRP  ,  ONLINE  ,  1ANOPT  ,  OPTION  ,  IFBACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

*  TYPE  DECLARATION 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,PURP0S(fo) 

LEVEL  2#  LU IN, LUOUT, LUDATA, LURUN, LUREPT, IFSVRP, ONLINE, I ANOPT 
LEVEL  2,  OPTION, IFBACH, USRNAM, PURPOS 

DATA  LU1N, LUOUT, LUDATA, LURUN, LUREPT  /  b,  6,  1,  2,  3  / 

*  SET  CONTROL  DEFAULTS. 

DATA  IFSVRP  /  0  / 

DATA  ONLINE  /0  / 

END 


PROGRAM  LISTING  (Continued) 
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o  o  o 


SUBROUTINE  AN AL 1 2 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*  [His  ROUTINE  CONTROLS  PROCESSING  OF  ALL  ANALYSIS  UPTIONS 
- FALCON  R&D  JLW  4/20/79 


*  SELECT  Aim  ANALYSIS  OPTION 


*  PROGRAM  CONTROL  VARIABLES 

COMMON  /  CONTRL  /  LU1N  ,  LUOUT  ,  LUOATA  ,  LURUN  ,  LUREPT 

COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH 
CUMMON  /  CONTRL  /  USKNAM  ,  PURPOS 


*  TYPE  DECLARATION 


INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,PURP0S(6) 

LEVEL  2,  LUIN, LUOUT, LUUAT A , LURUN, LUREPT, IFSVRP, ONLINE, IANOPT 
LEVEL  2  ,  OPTION, IFBACH, USRNAM, PURPOS 
C 

READ (LUIN, 9002)  IANOPT 
WRITE  (LUOUT,  9003)  IAimOPT 
C 

IF  (IANOPT  .60.  0)  GU  TU  9000 
1F(1Aim0PT.E0.3)CALL  SKPUAT 
IF  (IANOPT  .EO.  1)  CALL  KITTYtt 
IF  (IANOPT  .EQ.  2)  CALL  COMPLT 
IF  (IANOPT  .LT.  0)  CALL  PERROR(S) 

IF (IANOPT .GT.3ICALL  PERR0R(6) 

C 

9000  RETURN 


C 


9001  FORMAT (20X , 40HSELECT 

9002  FORMAT (12) 

9003  FORMAT(1H1,///,10X,1SHANALYSIS 


* 

* 

* 

* 

* 


/, 15X,23N0 
15X,28H1  = 
/,  15X,33H2 
/, 15X, 13H3 
//) 


ANALYSIS  OPTION  0,  1,  OR  2  N  U  R  . .  /) 


OPTION  IS, IS, 


=  NO  FURThER  ANALYSIS,/, 

USER  CONTROLLED  ANALYSIS, 

=  COMPLETE  DATA  RECORD  ANALYSIS, 
=  SKIP  DATA, 


END 


PROGRAM  LISTING  (Continued) 
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oooo  nnooonooo  ooo 


SUoKUUT  I  NE  COMPLY 


C - — - 

C 

C  *  SUBROUTINE  LUHPLT  PUNS  A  C0MPLE1E  SPECTRUM  OF 
C  TRANSFORMATIONS  ON  GUN-PROJECTILE  DYNAMICS  DATA 

C 

- - 

COMMON  /  CUMSUB  /  TIN1T  ,  HAST  ,  NODTS 

COMMON  /  SOATA  /  NSTPTS  ,  DTIME  ,  STHING(bOOO)  ,  STIME(5000) 

COMMON  /  DTECT  /  EP0CH1  ,  EPOCHS 

COMMON  /  REPDAT  /  REPSTR(5000)  ,  IFREP  ,  NUREP 

COMMON  /  CLASCM  /  ICLASS  ,  NBEGBL 

COMMON  /  REGII  /  MAXPT,  MINPT,  VALMAX,  VALMIN 

COMMON  /  CMPLT  /  TEND  1  ,  TENDS  ,  IUIV  ,  SUBTIM 

*  fTLENl  ,  TLENS  ,  T  LEN 
INTEGER  EPOCH  1  ,  EPOCHS 

LEVEL  S,  T IN  IT, TL AST, NODTS, NSTPTS, DTIME, STRING, ST IME, EPOCH  1, EPOCHS 

*  ,  ICLASS, NBEGBL, TEND1 , TENDS, IDIV, SUBTIM, NUREP 

*  ,  REPSTR, IFREP, TLEN1, TLENS, TLEN 

*  ,  MAXPT, MINPT, VALMAX, VALMIN 

C  *********************************************************************** 

c 

C  *  READ  IN  COMMON  DATA 

C  AND  DETERMINE  SATURATION  LEVEL 

C 

CALL  OVERLAY tlHb, 1, 01 
C 

C  *  CALCULATE  ThE  FIRST  AND  SECOND  DERIVATIVES  FOR 
C  UNSATURATION  PORTION  OF  DATASET 

C 

CALL  OVERLAY  (1HL, 7,0) 

C 

C  *  LOCATE  EPOCHS  1  AND  S  OR  EPOCH  1  IF  SATURATED  DATA  SET 
C 

CALL  OVERLAY  (1  HIM,  6, 0) 

C*********************************************************************** 

*  WRITE  PAGE  HEADING  FOR  COMPLY 
WRITE(6,9001 ) 

*  ANALYSIS  OVERLAYS  ARE: 

(C,S,0)  SUBSET 
( F , 5 , 0 )  THEFFT 
( D , 3 , 0 )  PAWS 
( P , 9 , 0 )  REGSAN 

*  CHOOSE  SUBSET  FUR  REGIME  ONE 

UNIT  =  ST  IME  Cl) 

TLAST  =  ST IME (EPOCH1 ) 

NODTS  =  1 

***************************  ******** 

*  CAN  REGIME  ONE  BE  SUBDIVIDED? 

TEND  1  =  SMALLEST  ALLOWABLE  TIME  FOR  SUBDIVISION  OF  REGIME  I 

C  TENDS  =  LARGEST  ALLOWABLE  TIME  FOR  SUBDIVISION  OF  REGIME  I 

PROGRAM  LISTING  (Continued) 
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ooo  ooo  o  on  o  o  o  o  ooo  ooo  o  o  o  o  o  ooo 


c 

T t i'll)  1  =  T  It'll  r  +  201  .  *  DTIME 
TEND2  =  TLAST  -  201.  *  DTIME 
IF  (IENU1  .GE.  T  E  N  D  2 )  GO  TU  1000 
3000  CUMT1NUE 

*  USER  DECIDES  IF  REGIME  ONE  NEEDS  10  BE  SUBDIVIDED 

WRITE(6,9002)  TINI1  ,  TLAST 
READ  (5,9003)  1 0 1 V 

*  I D I V  =  1,  ANALYSIS  CONTROL  RETURNS  TO  USER 

=  0,  REGIME  1  WILL  NOT  BE  SUBDIVIDED 
=  -1,  REGIME  I  WILL  BE  DIVIDED  bY  THE  USER 

IF  (1D1V  .ELI.  1)  GO  TO  9000 
IF  ( I D I V  .EQ.  0)  GO  TO  1000 
IF  (IDIV  .ED.  -1)  GO  TO  2000 
CALL  PERROR ( 7 1 ) 

2000  CONTINUE 

*  SUBTIME  IS  TIME  OF  SUBDIVISION  OF  REGIME  I 

READ  (5,9005)  SUBTIM 

*  CHECK  VALIDITY  OF  SUBTIME 

IF  (SUBTIM  .LE.  TINIT  .OR.  SUfaT IM  .GE.  TLAST) 

*  CALL  PERROR ( 7 3 ) 

IF  (SUBTIM  .GT.  TINIT  .AND.  SUbTIM  .LE.  TEND  1 ) 

*  SUBTIM  =  TENOl  +  DTIME 

IF  (SUBTIM  .GE.  TEND2  .AND.  SUbTIM  .LT.  TLAST) 

*  SUBTIM  =  TEND2  -  DTIME 

WRITE (b, 9004)  TEND  1  ,  TEND2  ,  SUbTIM 
***********  ********************* 

*  THREE  ANALYSES  OF  DA  1 A  -  TwO  SUBDIVISIONS  OF  REGIME  I 

NDREP  =3 

*  WRITE  HEADING  FOR  REGIME  I  PART  ONE 

IP  =  1 

TLAST  =  SUBTIM 

TLEN1  =  TLAST  -  TINIT 

WRlTE(fe,9006)  IP,  TINIT, TLAST, TLEN1 

*  SELECT  PART  ONE  OF  REGIME  ONE 

CALL  OVERLAY (1HC, 2, 0) 

*  ANALYZE  PART  ONE  OF  REGIME  ONE 

CALL  OVERLAY (1HF ,5, 0) 

CALL  OVERLAY (1HD, 3, 0) 

C 

C  *  WRITE  HEADING  FOR  REGIME  I  -  PART  TWU 
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c 


t. 


IP  =  2 

UNIT  =  SUttTIM 
T  L AST  =  STIME(EPUCHI) 

TLEN2  =  TLAST  -  UNIT 

WRITE  (6,  9006)  IP,  UNIT,  TLAST,  TLEN2 
C 

C  *  SELECT  PART  TWO  OF  REGIME  ONE 
C 

CALL  OVERLAY (1HC, 2, 0) 

C 

C  *  ANALYZE' PART  TWO  OF  REGIME  ONE 
C 

CALL  OVERLAY (1HF, 5, 0) 

CALL  OVERLAY (1HD,3,0) 

GO  TO  4000 
1000  CONTINUE 
C 

C  *  TWO  ANALYSES  OF  DATA 
C 

NDREP  =  2 
C 

C  *  CHOOSE  SUBSET  FOR  REGIME  ONE 
C 

CALL  OVERLAY (1HC,2,0) 

C 

C  *  WRITE  HEADING  FOR  REGIME  I 
C 

TLEN  =  TLAST  -  TINIT 
WRITE  16, 9007)  TINIT  ,  TLAST  ,  TLEN 
C 

C  *  RUN  THEFFT  AND  PAWS  ON  REGIME  ONE 
C 

CALL  OVERLAY (1HF, 5, 0) 

CALL  OVERLAY  ( 1HD, 3, 0) 

4000  CONTINUE 
C 

C  *  I  CL ASS  INDICATES  SATURATION  LEVEL: 


c 

I  CL ASS  = 

“1, 

HEAVY  SATURATION,  EPOCH  1  DETECTION 

c 

AND  REGIME  1  ANALYSIS. 

c 

ICLASS  = 

0, 

NO  SATURATION,  FULL  ANALYSIS 

c 

I  CL  ASS  = 

1  , 

LIMITED  SATURATION  (CLIPPING), 

c 

FULL  ANALYSIS. 

C 

IF  (1CLASS  .EG.  -1)  GO  TO  9000 

C***  *** ***************************************************************** 
C 

C  *  ANALYSIS  OF  REGIME  II 
C 

C  *  CHOOSE  A  SUBSET 
C 

TINIT  =  ST IME (EPOCH1 ) 

TLAST  s  STIME (EP0CH2) 

NODTS  =  1 

CALL  OVERLAY (1HC,2,0) 

C 

C  *  WRITE  HEADING  FOR  REGIME  II 
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c 

TLEN  =  TLAST  -  UNIT 
WRIT E (b, 9008)  UNIT,  TL  AST  ,  TLfcN 

*  SPECIAL  REGIME  II  ANALYSIS 
CALL  OVERLAY (IMP, 9,0) 

WRITE  (6,9010)  MAXPT,  VALMAX,  MINPT,  VALMIN 

*  POWER  SPF-CTkUM 


CALL  OVERLAY (1HD,3»0) 


;  *  * 


*  CHOUSE  SUBSET  FOR  REGIME  I H K E E 

UNIT  =  ST  I  ME  (EP0CH2) 

TLAST  =  STIME(NSTPTS) 

NOUTS  =  1 

CALL  OVERLAY (1HC, 2, 0) 

*  WRITE  HEADING  FOR  REGIME  III 

TLEN  =  TLAST  -  UNIT 

WRITE  (6,  9009)  UNIT,  TLAST,  TLEN 

*  RUN  THEFFT  AND  PAWS  ON  REGIME  THREE 


CALL  OVERLAY (1HF,5,0) 
CALL  OVERLAY ( 1HD, 3, 0) 


I******** 


9000  CONTINUE 

90  01  FORMAT  (1H1,///,20X,54(1H*),/,20X, 1H*,52X, 1H*,/,20X, 

*  lH*,9X,34hANALYSIS  OF  REGIMES  I,  II,  AND  III,9X, 

*  1H*,/,20X,1H*,S2X,1H*,/,20X,54(1H*)) 

9002  FORMAT (IX,////, 10X,38HREG1ME  ONE  IS  DEFINED  TO  START  AT  TIME, 

*  FS.a,2X, 15HAND  END  AT  T IME , F8 . 4 , 2X , 1 H . , / ,  1  0 X , 

*  bOHIF  IT  HAS  TWO  DISTINCT  SECTIONS  TU  IT,  1  HEN  IT  IS  , 

*  31HBEST  TO  DIVIDE  IT  IN  TwO  PARTS.) 

9003  FORMAT (12) 


9004  FORMAT (IX,// / , 10X,44HTHE  OPTION  TO  SUBDIVIDE  REGIME  ONE  HAS  BEEN  , 

*  8HCH00SEN. , /, 10X, 

*  4SHTHE  POINT  OF  SUBDIVISION  MAY  BE  BETWEEN  TIME, 

*  F 8 . 4 , 9H  AND  T I ME , F 8 . 4 , / , 1 0 X , 

*  3bH  THE  POINT  OF  SUBDIVISION  IS  AT  T I  ME , F8 . 4 , 2X , 1 H . ) 

900b  FORMAT (F8.4) 

9006  FORMAT (1H1,////,20X,30(1H*)»/»20X,1H*,28X,1H*,/,20X, 

*  1H*,5X, 1 7HREG1ME  I  AN AL Y S I S , 6 X , 1 H* , / , 

*  20X, 1H*,28X,  1  H*  , / ,20X,30(1H*) ,///,20X,4HPART ,12,/, 

*  30X,F8.4,2X, lbHMILLISECONDS  T0,F8.4,2X, 

*  12HMILLISEC0NDS, /,30X,F8. 4, 2X, 1 7HM ILL  I  SECONDS  LONG) 

9007  FORMAT ( 1 H 1 ,////, 20X , 30 ( 1 H*), /, 20X, 1 H* , 28X , 1 H* , / , 

*  20X, 1H*,5X, 17HREGIME  I  AN AL YS I S , bX , 1 H* , / , 

*  20 X,  1H*,28X, 1H*,/,20X,30(1H*) ,///,30X, 

*  F8.4,2X, lbHMILLISECONDS  T0,F8.4,2X, 

*  12HMILLISECONDS,/,30X,F8.4,2X,  17HM1LHSEC0NDS  LONG) 
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90  0  a  FURMAT(1H1,////,2UX,30(1H*),/»<J0X,1H*,26X,1H*,/,I?0X,1H*,SX, 

*  lHhwEGlME  II  ANALYSIS, bX, 1H* , /,20X, 1H*,26X, 1H*, /,20X, 

*  30llH*),///,30X,Ftt.4,2X, lbHMlLLl SECONDS  lO,Ftt.4,2X, 

*  12HMILL1SECUNUS, /,30X,Fb.4,2X» 17HMILLISECUNUS  LONG) 

90  0  9  FORMAT  (1H1,////,20X,3011H*),/,20X,  1H*,28X,1H*,/,20X,1H*,5X, 

*  19HKEGIME  III  ANAL YSIS# OX , 1 H* , /, 20X, 1H* , 2SX # 1H*, / , 20X# 

*  30 ( 1H*) , ///, 30X,F6.4, 2X, 1 bhKl LL 1  SECONDS  T0,F8.4,2X, 

*  12HMlLLISEC0NDS,/,30X,Ftt.4,2X, 1 7 HM ILL  1  SECONDS  LONG) 
9010  FORMAT (IX,////, 1 0 X , 25HM A  X  I  MUM  ACHIEVED  AT  POINT, 

*  1 1  0 , 2  X  ,  lOHwlTH  VALUE, F10. 0, /, 10X, 

*  2bH MINIMUM  ACHIEVED  AT  POINT , 1 1 0 , 2 X , 

*  IOHnITH  VALUE,  F10.4,///) 

RETURN 

END 
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SUBROUTINE  1AINP1 


*  THIS  ROUTINE  PROMPTS  INTERACTIVE  INPUT  FROM  THE  USER. 
- FALCON  R&U  JLW  6/5/79 


*  PROGRAM  CUNTROL  VARIABLES 

COMMON  /  CUNTKL  /  LUIN  ,  LUOUT  ,  LUOATA  ,  LURUN  ,  LUREPT 

COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFbACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

*  TYPE  OECLARATION 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM  ( 2 )  ,Pl)RP0S(b) 

LEVEL  2,  LUIN, LUOUT, LURUN, LUREPT, LUO AT A, IFSVRP, ONLINE, IANOPT 
LEVEL  2,  OPTION, IFBACH, USRNAM, PURPOS 

*  PROMPT  USER  NAME  INPUT 

WRITE(LU OUT, 9001) 

C 

REAOILUIN,*)  USRNAM 
WRITE (LUOUT , 9005)  USRNAM 
C 

C  *  PROMPT  TEXT  TO  DOCUMENT  PURPOSE  OF  THE  RUN. 

C 

WRITE (LUOUT ,9002) 

C 

CALL  ROPURP 
C 

9000  RETURN 
C 

9001  FORMAT (IX, ////, 10X, 1 1HUSER 'S  NAME) 

9002  FORMAT (IX,///, 10X, 19HPURPOSE  OF  RUN) 

9005  FORMA  T (15X,2A10) 

C 

END 
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SOHROU  T  I  NE  InCIRL 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*  IHIS  ROUTINE  CONTROLS  INPUT  RtADING  AND  INTERPRETATION  OF  USER 

*  RUN  CONTROLS. 

- — - - - FALCON  R SO  JLW  4/20/7R 


*  PROGRAM  CONTROL  VARIABLES 

COMMON  /  CONTRL  /  LUIN  ,  LUOIJT  ,  LUOATA  ,  LURUN  ,  LUREPT 

COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

*  TYPE  DECLARATION 


INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,PURPOS(b) 

LEVEL  2, LUIN,  LUOUT,  LUD AT A,  LURUN, LUREPT, IFSVRP, ONLINE, IANOPT 
LEVEL  2,  OPTION, IFBACH, USRNAM, PURPOS 

*  DETERMINE  PROCESSING  MODE 


*  BATCH  MODE  REQUIRED 

IFBACH  =  0 
ONLINE  =  0 
CALL  IAIiMPT 
C 

END 
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SUBROUTINE  INFILE 


THIS  ROUTINE  PREPARES  ALL  NECESSARY  FILES  FOR  THE  RUN. 

- - - FALCON  K NO  JLW  4/20/7R 


PROGRAM 

CONTROL  VARIABLES 

COMMON 

/ 

CUNIRL 

/ 

LUIN  , 

LUOUT  ,  LUOATA  ,  LURUN 

,  LUREPT 

common 

/ 

CONTRL 

/ 

IFSVRP 

,  ONLINE  ,  IANOPT  ,  OPT 

ION  ,  IFBACH 

COMMON 

/ 

CONTRL 

/ 

USRNAM 

,  PIJRPOS 

*  TYPE  OECLAKATION 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM 12)  ,PURP0S(6) 

LEVEL  8,  LUlN, LUOUT , LUO AT  A, LURUN, LUREP1 , IFS VRP , ONL I NE , 1 ANOPT 
LEVEL  a,  OPTION, IFbACH, USRNAM, PUkPOS 


REAINO  LUO A 1  A 


*  RESERVE  PF  OEVICE 


FOR  REPORT  FILE 


9000  RETURN 
ENO 
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SUBROUTINE  K1TTY8 

C - - - - - - - 

C 

C  *  SUBKOUT ] NE  RITTY8  CONTROLS  THt  ANALYSIS  OF  GUN  DATA, 

C  BY  SELECTING  VARIOUS  OPTIUNS 

C 

C - FALCON  R&D  MEK  6-19-79 

COMMON  /  CONTRL  /  LUIN  ,  LUOUI  ,  LUOATA  ,  LURUN  ,  LUREPT 
COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACM 
CUMMON  /  CONTRL  /  USRNAM  ,  PURPOS 
COMMON  /COMOPT/  SUBOPT 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,  PURPOS(b) 

INTEGER  SUBOPT 

LEVEL  2,  LUlN,LUOUT, LUOATA, LURUN, LU REP  I , ONLINE, IFSVRP, IANOPT 
LEVEL  2, OPT  ION, USRNAM, PURPOS, SUbOPT 
LEVEL  2,  IFBACH 
C 

C  *  Rt AD  IN  COMMON  SOA1A 
C 

CALL  OVERLAY  (1F)B,  1,0) 

1000  CONTINUE 

READ  (LUIN  ,9002)  SUBOPT 
3000  CONTINUE 

WRITE  (6,9003)  SUBOPT 
IF  (SUbOPT  .ED.  0)  GO  10  9000 
IF  (SUBOPT  .EG).  11)  GO  TO  3500 
IF  (SUBOPT  .EQ.  12)  GO  TO  3750 

IF  (SUBOPT  .LI.  1  .UK.  SUBOPT  .GT.  5)  CALL  PERR0K(79) 

Gu  TO  (9000,  5000,  6000,  7000,  7500),  SUBOPT 
GO  TO  9000 

C*ft*A**A*******A****ft**ft*ft****** 

3500  CONTINUE 
C 

C  *  SELECT  ANALYSIS  0P1 ION  OERIV  --  FIRST  AND 
C  SECOND  DERIVATIVES  OF  ENTIRE  DATA  SET. 

C 

CALL  OVERLAY ( 1HL, 7, 0) 

GO  TO  8000 

C******************************** 
3750  CONTINUE 
C 

C  *  SELECT  ANALYSIS  OPTION  DETECT  --  FIND  FIRST  AND 
C  SECOND  EPOCHS  IN  DATA  SET. 

C 

CALL  OVERLAY ( 1HN, 8, 0) 

GO  TO  8000 

C******************************** 

9000  CONTINUE 
C 

C  *  SELECT  ANALYSIS  OPTION  SUBSET  --  CHOOSE  SUBSET  OF  DATA. 

C 

CALL  OVERLAY ( 1HC, 2, 0) 

GO  TO  8000 

C********************************* 
5000  CONTINUE 
C 

C  *  SELECT  ANALYSIS  OPTION  PAWS  --  POWER  SPECTRUM  ANALYSIS 
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CALL  OVERLAY (1HD, 3, 0) 

GU  10  8000 

C  *  *  *********************  ************ 

6000  CONTINUE 

*  SELECT  ANALYSIS  OPTION  GIRAFfc  --FORWARD  BACKWARD  LEAST  SQUARES  FITS 

NOT  CURRENTLY  AVAILABLE 

CALL  OVERLAY (1HE,4,0) 

GO  TO  8000 

*********************************** 

7000  CONTINUE 

*  SELECT  ANALYSIS  OPTION  THEFFT  --  CREATES  «TH  DEGREE 

POLYNOMIAL,  PERIODIC  FUNCTION  OF  FIRST  DEGREE 
FREQUENCIES  TO  FIT  DATA  SUBSET. 

CALL  OVERLAY (1HF,5,0) 

GO  TO  8000 

*********************************** 

7500  CONTINUE 

*  SELECT  ANALYSIS  OPTION  SMOOTH  —  11-POINT  LOCAL 

AVERAGE  OR  11-POINT  LOCAL  MEDIAN. 

CALL  OVERLAY(1HG,6,0) 

GO  10  8000 

C*********************************** 

8000  CONTINUE 

GO  TO  1000 
9000  RETURN 

9003  FORMAT  (  1H1 ,////, 10X, 19HANALYSIS  OPTION  IS  ,I2,/,15X, 

*  1 1 H 1  =  SUBSET,/, 15X, 

*  9H2  s  PAWS, /, 15X, 1 1H3  =  GIKAFE,/, 15X, 1 1  H#  =  THEFFT, 

*  / , 1 5X » 1 1 H5  =  SMOOTH,/, 15X,10H11  =  UER1V,/,15X, 

*  1 1  H12  =  DETECT) 

9002  FORMAT  (12) 

END 
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SUBROUTINE  PERRUR  UERknO) 


C 

c 

c 

c 


*  ERROR  REPORTING  ROUTINE  -  SEE  USER  DIAGNOSTICS  FOR  MORE  INFORMATION. 
- -FALCON  k&D  JLW  4/20/79 


*  PROGRAM  CONTROL  VARIABLES 

COMMON  /  CONTRL  /  LUIN  ,  LUOUT  ,  LUO A  T A  ,  LURUN  ,  LUREPT 

COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  1AN0PT  ,  OPTION  ,  IFBACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

*  TYPE  DECLARATION 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM 12)  ,PURPOS(b) 

LEVEL  2,  LU IN, LUOUT, LURUN, LUD AT  A , LUREP T , I F S VRP , ONL I NE , I ANOPT 
LEVEL  2,  OPTION, IFBACH, USRNAM, PURPOS 

WRITE(LUUUT,9001)  IERRNO 

*  SPECIAL  ERROR  PROCESSING  GOES  HERE. 

9000  STOP  "  ABNORMAL  TERMINATION  ***  SEE  OUTPUT  AND  USER  DIAGNOSTICS**" 
C 

9001  FORMAT  (IX,  19HPR0GRAM  ERROR  TYPE  ,13,  19H  HAS  BEEN  DETECTED  ,/, 

1  1  X , 4 1 HSEE  USER  DIAGNOSTICS  FOR  MORE  I NF ORMAT I  ON , 1 0 ( 1 H* ) ) 

C 

END 
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SDURUUT  I  ME  ROPURP 


C 

*  THIS  ROUTINE  ACCEPTS  USER  INPUT  TO  DOCUMENT  THE  PURPOSE  OF  THE  RUN. 

*  THE  TEXT  IS  WRITTEN  TO  A  RUN  INFU  FILE  TAPES  (LURUN  =  2) 

- - - FALCON  R&D  JLW  b/5/79 


*  PROGRAM  CONTROL  VARIABLES 

COMMON  /  CONTRL  /  LUIN  ,  LUUUT  ,  LUDATA  ,  LURUN  ,  LUREPT 

COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFHACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

*  TYPE  DECLARATION 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM  1 2 )  ,PURPUS(b) 

LEVEL  2,  LU I N,LUOUT, LUDATA, LURUN, LUREPT , I FS VRP , ONL 1 NE , IANOPT 
LEVEL  S,  OPTION, IFBACH, USRNAM, PURPOS, 1LINES, NL1NES 
COMMON  /  REPCUM  /  ILINES  ,  NL1NES 

REAU(LUlN,9001)  NL1NES 
DO  3000  ILINES=1 ,N LINES 
1000  READ(LU1N,9003)  (PURPOS (I),  1  =  1,6) 

WRITE (LUOUT, 9002)  ( PURPUS II ) , 1 = 1 , 6 ) 

WRITE (LURUN, 9002)  (PORPOS(I) , 1=1,6) 

3000  CONTINUE 

*  END  OF  INPUT  TEXT,  REWIND  FILE  AND  EXIT. 

7000  REWIND  LURUN 
C 

9000  RETURN 
C 

9001  FORMAT (15) 

C 

9002  FORMAT (15X,6A10) 

9003  FORMAT (6A10) 

C 

END 
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SUBklJUl  1  Nt  SKPDAT 

COMMON/CONTkL/LUlN,LUOur,LUDAlA,LUkUN,LURtPT 
COMMON/CONT KL/ IF SVRP, ONLINE, 1ANOPT , OPTION »IFbACH 
COMMON /CON  TKL/IJSRN  AM,  PURPOS 
INTEGER  ONLINE,  OPTION,  USRNAM(2),PURP0S(6) 

LEVEL  2, L UTN,LUOUT,LUDATA,LURUN,LUkEPT, ONLINE, IFSVRP,IANUPT 
LEVEL  2, OPTION, USRN AM, PURPOS, IFBACH 

READ  DATA  ONLY  -  DON'T  PROCESS 
ILOOP=  NO.  OF  CHANNELS  TO  BE  SKIPPED 

READ(LUIN,9001) ILOOP 
DO  1000  IL=1, ILOOP 
CALL  OVERLAY ( 1HB, 1 , 0) 

1000  CONTINUE 

READILUIN, 9002) I ANOPT 

9000  RETURN 

9001  FORMAT (15) 

9002  FORMAT (12) 

END 
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SUBKUUT 1NE  START 


L - - - - - — - - 

C 

C  *  ThIS  ROUTINE:  HANDLE: S  SIAkT  UP  PROCESSING. 

C 

C - FALCON  R&U  JLh 

C 

C 

C  *  PREPARE  FILES  FOR  THE  RUN. 

C 

CALL  INFILE 

*  READ  AND  INTERPRET  USER  CONTROLS  FUR  THIS  RUN 
CALL  1NCTRL 

*  OTHER  REQUIRED  INITIALIZATION  FOR  THIS  RUN. 

RETURN 
END 
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OVERLAY  IB, 1  ,  U) 
PROGRAM  DAT  AIN 


*  SlIBROUT INE  UATAlN  READS  MUZ/LEUE  V ICE  DATA  FROM  TAPE1 

- - - - - FALCON  R&D  NLA  b-7-/9 

COMMON  /  SDaTA  /  NSTPTS  ,  DTlMt  ,  STRlNG(bOOO)  ,  STlMt(5000) 

COMMON  /  LABEL  /  KTB124)  ,  TST APT  ,  I 

COMMON  /  REPDAT  /  REPSTR  IbOOO)  ,  IFREP  ,  NDREP 

LEVEL  2,  REPSTR,  IFREP,  NDREP 

LEVEL  2,  NSTPTS,  DTIME,  STRING,  STIME,  KTB,  I,  T  ST  ART 

*  READ  AND  WRITE  LABEL 

READ  III  KTB 
WRITE  lb,  9001)  KTB 

*  READ  AND  WRITE  NSTPTS  =  NUMBER  OF  DATA  POINTS 

TSTART  =  START  TIME  OF  ARRAY 
DTIME  =  TIME  INCREMENT  OF  ARRAY 

READ  (1)  NSTPTS 

WRITE  lb,  9002)  NSTPTS 

READ  (1)  TSTART,  DTIME 

WRITE  lb,  9003)  TSTART,  DTIME 

*  READ  DATA  IN  ARRAY  STRING 

READ  11)  (STRINGll),  1=1, NSTPTS) 

*  SET  UP  TIME  ARRAY  STIME 

STIMEll)  =  TSTART 
DO  1000  1=2, NSTPTS 
STIME II)  =  STIMEll  -  1)  +  DTIME 
1000  CONTINUE 

*  PLACE  DATA  STRING  INTO  ARRAY  REPSTR.  IFREP  INDICATES  THE 

REGIMES  REPLACED  BY  THEFFT  INTO  REPSTR. 


IFREP  =  0 

DO  2000  I  =  1, NSTPTS 
REPSTRII)  =  STRINGll) 
2000  CONTINUE 


C 

CALL  CLASIF 


C 


9000  RETURN 

9001  FORMAT UHl,///,20X,b4(lH*),/,20X,lH*,52X,lH*,/,20X,lH*, 

*  14X,24HMUZZLEDEVICE  DATA  RECORD , 1 0 X , 1 H* ,/, 20 X , 


★ 

★ 

9002  FORMAT 

9003  FORMAT 

* 


lH*,52X,lH*,/,20X,lh*,2X,24I2,2X,lH*,/,20X,lH*, 

b2X,lH*,/,20X,54(lH*),//) 

110X,  2 1 HNUMBER  OF  DATA  POINTS,  110,/) 

110X,  12HIN1T IAL  TIME,  F10.4,  bX,  14HTIME  INCREMENT, 
10X,  22HITIME  IN  MILLISECONDS)) 


END 


F10.4, 
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SUBROUT 1 NE  CLASIF 


*  SUBROUTINE  CLASIF  CLASSIFIES  THE  DATA  SET 
ACCORDING  TO  THE  LEVEL  OF  S A TUk A T 1  ON . 

C 

C - - - 

LEVEL  2,  NLIMBL,  NBLACK,  NENOBL,  <JS  TAR,  NSTPTS, STRING, NBEGBL,  I  CL  ASS,  1 
LEVEL  2,DTIME,ST IME 

COMMON  /  SO A T A  /  NSTPTS  ,  DT IME  ,  STRING(SOOO)  ,  STIMEC5000) 
COMMON  /  CMCLAS  /  NLIMBL,  NBLACK,  NENOBL ,  QSTAR,  I 
COMMON  /  CLASCM  /  ICLASS,  NBEGBL 
C 

C  *  DETERMINE  MAXIMUM  VALUE  OF  DATA  RECOkU  rQSTAk 

C  FIRST  OCCURRENCE  =NBEGBL 

C  NUMBER  OF  OCCURRENCES  sNBLACK 

C 

C  NBLACK  =  0  NO  SATURATION 
C  ICLASS  =  0 

C 

C  NBLACK  <  10  LIMITED  SATURATION 
C  ICLASS  =1 

C 

C  NfaLACK  >  10  HEAVY  SATURATION 
C  ICLASS  =  -1 

C 

NLIMBL  =  10 
NBLACK  =  0 
NENOBL  =  0 
NBEGBL  =  0 
QSTAR  =  -99 999 . 

DO  5000  1  =  1, NSTPTS 

IF  (STRING(I)  -  QSTAR)  1000,  2000,  3000 
1000  CONTINUE 

GO  TO  9000 
2000  CONTINUE 

NENOBL  =  I 
GO  TO  9000 
3000  CONTINUE 

QSTAR  =  STRING (I) 

NBEGBL  =  I 
9000  CONTINUE 
5000  CONTINUE 

NBLACK  =  NENDBL  -  NBEGBL 
ICLASS  =  0 

IF  (NBLACK  .LE.  0)  NbLACK  =  0 
IF  (NBLACK  .LE.  0)  GO  TO  bOOO 
ICLASS  =  -1 

IF  (NBLACK  .LT.  NLIMBL)  ICLASS  =  1 
6000  CONTINUE 

WRITE  (6,9001)  ICLASS,  NBLACK 
IF  (ICLASS  ,EQ.  0)  WRITE (b, 9002) 

IF  (ICLASS  .EG).  1)  WRI  TE(b,9003) 

IF  (ICLASS  .EG.-l)  WRITE(6,9009) 

9000  RETURN 

9001  FORMAT (IX, /////, 20X,27HT1ME  HISTORY  CLASSIFICATION,//, 

*  10X,24HCLASSIFICATION  NUMBER  IS, 15,/, 
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*  10X, 35HNUMBER  OF  POINTS  IN  SIGNAL  BL ALKOUT , I  1 0 ) 

900d  FORMAT  (1  X,  /  ,  1  OX  ,  cf9HC0MPLt  T  E  ANALYSIS  OF  DATA  SET, 

*  13H  IS  PERMIT ftO,//) 

9003  FORMAT ( 1 X ,/, 1  OX, 31 HLlMl TED  SAlUkATlON  UF  DATA  SET,, 

*  3 1 H  COMPLETE  ANALYSIS  IS  PERMITTED,//) 

9009  FORMAT (IX,/, 10X,c?9HHEAVY  SATURATION  OF  DATA  SET,, 

*  55H  ONLY  EPOCH  1  DETECTION  ANU  REGIME  1  ANALYSIS  PERMITTED, 

*  //) 

END 
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PROGRAM  SUbSET 


c 

C  *  SUBROUTINE  SUBSET  TAKES  A  SUBSET  OF  THE  DATA  ARRAY  STRING  FOR 
C  FURTHER  ANALYSIS 

£ _ FALCON  R&D  NEA  b-5-79 

COMMON  /  SOATA  /  NSTPTS  ,  UTIME  ,  STRING(5000)  *  STIME(SOOO) 
COMMON  /  SSOATA  /  NSSPTS  r  SDTIME  ,  SUBSTR(5000)  ,  SSTIME(5000) 
LEVEL  2,  NSTPTS,  DTIME,  STRING,  STIME,  NSSPTS,  SDTIME, 

*  SUBSTK ,  SSTIME 

LEVEL  2,  T 1NIT , TLAST ,NOOTS, IS, IE, INO, I, LIMP 
COMMON  /  COMSET  /  I S , 1 E , I ND , 1 , L 1 MP 
COMMON  /  COMSUB  /  TINIT  ,  TLAST  ,  N0DT5 
C 

C  *  READ  parameters 

TINIT  =  INITIAL  TIME  DESIRED  FOR  SUB-STRING 
TLAST  =  FINAL  TIME  DESIREO  FOR  SUB-STRING 

GOUTS  =  NUMBER  OF  DTIME  INTERVALS  BETWEEN  SUB-STRING  POINTS 
CALL  SUbPAR  (TINIT,  TLAST,  NOOTS) 

*  INITIAL  AND  FINAL  INDICES  IN  STRING 

IS  =  (TINIT  -  STIME  (1JT/DTIME  +  1 

IE  =  (TLAST  -  ST1ME(1) )/DTIME  +  1 

*  CHECK  INDICES  OF  STRING  ARRAY  TO  MATCH  TIMES 

IF  (IS  .LE.  0)  CALL  PERR0R(53) 

IF  (IE  .LE.  1ST  CALL  PERRQK(5«) 

IF  ( ABS (ST IME ( I S)  -  TINIT)  .GT.  DTIME)  CALL  PERROR(feO) 

IF  ( ABS (ST IME ( IE)  -  TLAST)  .GT.  DTIME)  CALL  PERROR(bl) 

IF  (TINIT  .EO.  STIME(IS+1))  IS  =  IS  +  1 
IF  (TLAST  .EO.  ST  I  ME ( 1 E  + 1 ) )  IE  =  IE  +  1 
C 

C  *  CREATE  SUBSTR  AND  SSTIME  ARRAYS  FROM  STRING  AND  STIME 
C 

I  =  0 

DU  1000  1N()=IS,  IE,  NODTS 
1  =  1  +  1 

SUBSTR ( I )  =  STRING(IND) 

SSTIME(I)  =  STIME(IND) 

1000  CONTINUE 
C 

C  *  SET  NSSPTS  TO  NUMBER  OF  POINTS  IN  SUBSTR  AND  SSTIME 
C 

NSSPTS  =  1 

IF  (NSSPTS  .LE.  0)  CALL  PERR0R(55) 

C 

C  *  SET  SDTIME  TO  TIME  INTERVAL  IN  SUBSTR  AND  SSTIME 
C 

SDTIME  =  DT1ME+FL0AT (NODTS) 

C 

9000  RETURN 
END 
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SUBROUT  I  NE  SUBPAR  (UNIT,  TLAST ,  NODTS) 


*  SUBROUTINE  SUBPAR  READS  1NI1IAL  AND  FINAL  TIMES  WITH  POINT  INIERVAL 

FOR  CHOOSING  SUBSET  OF  SAMPLE  DATA 

- - - - - - FALCON  K&D  NEA  6-5-7* 

COMMON  /  CONTKL  /  LUIN  ,  LUOUT  ,  LUOATA  ,  LUKUN  ,  LUREPT 
COMMON  /  CONTRL  /  1FSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPUS 
INTEGER  ONLINE  ,  OPTION  ,  USRNAM (2)  ,  PURPOS(b) 

LEVEL  2,  LUIN, LUOUT, LUD AT  A , LURUN , LUREPT, ONLINE, 1FSVRP, IANOPT 
LEVEL  2,  OPTION, USRNAM, PURPOS , T IN  1 T , TLAST, NODTS 
LEVEL  2, IFBACH 

*  READ  TINIT  =  INITIAL  TIME 

TLAST  =  FINAL  TIME 

NODTS  =  NUMBER  OF  UTlME'S  BETwEEN  SUBSTRING  POINTS 


IANOPT  =  2  ANALYSIS  OPTION  IS  COMPLT 

UNIT,  TLAST, NODTS  ARE  SPECIFIED  PRIOR  TU  SR  CALL 

IF  (IANOPT  .EU.  2)  GO  TO  1000 
READ  (LUIN,  9001)  UNIT,  TLAST,  NODTS 
1000  CONTINUE 

IF  (TLAST  .LE.  UNIT)  CALL  PERR0R(51) 

IF  (NODTS  .LE.  0)  CALL  PERROR(52) 

IF  (IANOPT  .EU.  2)  GO  TO  2000 
WRITE  (LUOUT,  9002)  UNIT,  TLAST,  NUDTS 
2000  CONTINUE 

9000  RETURN 

9001  FURMAT  (2F10.4,  110) 

9002  FORMAT (IX,/, 2BH  DATA  SUBSET  BEGINS  AT  TIME  ,F10.A, 

*  28H  MILLISECONDS,  ENDS  AT  TIME  ,  F10.4, 

*  37H  MILLISECONDS,  AND  TAKES  POINTS  EVERY,  lb,  7H  DT1MES,/) 
END 
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c- 

c 

c 

c 

c 

c 

c 

c- 


OVERLAY  CD, 3, 0) 
PH OGHAM  PAWS 


*  SUBROUTINES  PAWS  CALCULATES  THE  PuwEK  SPECTRUM  FOR  ARRAY  SUBSTR 
USING  SYSTEM  SUBROUTINE  FTFREU 

CALLED  BY  COMPLY  (1AN0PT=2)  OR  K1TTY8  (IANOPT=l) 

PAWS  =  POWER  ANALYSIS  OF  WAVEFORM  SPECTRA 


C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 


COMMON  /  CONTRL  / 


COMMON 

COMMON 


CUMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/  S50ATA 
/  PAWCOM 


/FREOPS/ 

/FNEUS/ 

/  SCRIB1 
/  LAGCOM 
/  PAW VAR 


- FALCON  R&U  NEA  b-b-/' 

LUIN  ,  LUOUT  ,  LUUATA  ,  LURUN  ,  LUREPT, 

IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH, 
USRNAM  ,  PURPOS 

/  NSSPTS  ,  SOUME  ,  SUBSTR  15000)  ,  SST IME  ( 500  0 ) 

/  XYMV(b),  AC  V  ( 1 b02) ,  FREU(BOl),  PS(lb02), 

XCOV ( lb02) r  XSPECT l lb02) ,  AMPHAS ( 1 b02) , 

XFER  C lb02) ,  COHER (SOI ) 

A ( 8 1 0  J  ,  B(810)  ,  R(lb20)  ,  NROW 

FREQENC20)  ,  NFR 

/  XARR(IOO)  ,  YARRC100)  ,  UXI  ,  NDPTS  r  NPP1 

/  M 

/  XC1-500)  ,  INOlb)  ,  X I  NO  12  J  ,  PLOTOP  ,  LAGOPT, 

I,  NFREO 


INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,  PURPOS(b)  ,  PLOTOP 

LEVEL  2,  NSSPTS,  SOT  I ME ,  SUBSTR,  SS1IME,  A,  B,  R,  FRE9FN,  NFR 

*  ,  LUIN, LUOUT, LUDA1A, LURUN, LUREPT, IFSVRP, ONLINE, IANOPT 

*  ,  OPTION, USRNAM, PURPOS, M,NROW 

*  ,  IFBACH, XYMV, ACV,FREO, PS, XCOV, XSPECT, AMPHAS, XFER, COHER 

*  ,  X, INO, XIND, PLOTOP  ,  LAGOPT , I , NFREO 

*  WRITE  PAGE  HEADING  FOR  PAWS 
WRITE (b, 9005) 

PLOTOP  =  PLOTTING  OPTION  FOR  POWER  SPECTRUM 
=  0  OFF 
=  1  ON 

LAGOPT  =  FREUUENCY  MULTIPLE  OPTION  FOR  POWER  SPECTRUM 

=  0  MINIMUM  ALLOWABLE  FREUUENCY  MULTIPLE  1.25,. 5,1  OR  2) 
=  1  FREUUENCY  MULTIPLE  IS  1.0 

READC5,9003)  PLOTOP  ,  LAGOPT 
wR!TE(b,9004)  PLOTOP  ,  LAGOPT 


PARAMETERS 

M  =  NUMBER 


OF  LAG  TIMES 


500 


M  = 

2.  *  (1 

.  /  SOT IME) 

IF 

(LAGOPT 

.EU.  1) 

GO  TO 

5 

00 

IF 

(M  .GT. 

NSSPTS) 

M  =  M 

/ 

2 

IF 

(M  .GT. 

NSSPTS) 

M  =  M 

/ 

2 

IF 

CM  .GT. 

NSSPTS) 

M  =  M 

/ 

2 

GO 

TO  7  50 

CONTINUE 
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M  =  M  /  4 

IF  (M  .  G  T  .  NSSPTS)  M  =  M  /  2 
IF  ( M  .  b  T  .  NSSPTS)  M  =  M  /  2 
C 

C  *  INITIALIZE  F  T  FREQ  ARRAYS 
C 

750  CONTINUE 

DO  1000  1=1, NSSPTS 

x(i)  =  substru) 

1000  CONTINUE 

I  NO  ( 1 )  =  0 

IND(2)  =  NSSPTS 
1ND(3)  =  0 

IND(4)  =  M 
IN0C5)  =  0 
I  NO (6)  =0 

XINO(l)  =  SDTIME 
XIND(2)  =  0 
C 

C  *  CHECK  INPUT  TO  SUBROUTINE  FTFREU 
C 

IF  ( INI)  ( 2)  .LE.  3)  CALL  PERNOR  157  ) 

IF  ( I NO ( 0 )  .LE.  2)  CALL  PERKOR  (56) 

IF  (  INI)  (4)  .GE.  I  NO  C  2 )  )  CALL  PERRON  (59) 

*  00  SPECTRUM  ANALYSIS 

CALL  FTFREU  (X,  IND,  XINO,  XYMV,  ACV,  FREQ,  PS, 

*  XCOV,  XSPECT,  AMPHAS,  XFER ,  COHER,  1ER) 

IF  ( 1ER  .NE.  0)  WRITE  (LUOUT,  9001)  IER 

*  PRINT  OUTPUT  FRUM  POWER  SPECTRUM 

CALL  PAwOUT 

*********************************** 

*  PUT  TOP  10  MAJOR  FREQUENCIES  (FREQUENCIES  WITH  LARGEST 

POWER  SPECTRUM)  INTO  FREQEN. 

NFR  =  10 
NROW  =  M  +  1 
DO  3000  I  =  1 , NROW 
A(I)  =  FREQ(I) 

B ( I )  =  PS(I) 

3000  CONTINUE 
C 

C  *  SORT  FREQUENCIES  INTO  R  BY  DECREASING  POWER 
C 

CALL  ROWSKT (A, B,R,NROw) 

C 

C  *  PUT  10  LARGEST  (FIRST  10)  INTO  FREQEN 
C 

DU  4000  I  =  1 , NFR 
FREQEN ( I )  =  R ( 1 ) 

4000  CONTINUE 
C 

C  *  WRITE  10  PREDOMINANT  FREQUENCIES 
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WRITE  (6, 9002)  (FREQEN(I)  ,  1  =  1  » NFk  ) 

************************  *********** 

IF  ( PLOTuP  .EG).  0)  GO  TO  7000 

*  PLOT  LOG(POwER)  VS  FREQUENCIES  (FIRST  25) 

NFREO  =  H  +  1 
DO  bOOO  I  =  1, NFREO 

IF  (PSCI)  .LT.  1)  PS(I)  =  1. 

PS(I)  =  ALOGIOIPS(I)) 

5000  CONTINUE 

0X1  =  FREU(2)  -  FREQ ( 1 ) 

NOPTS  =  MINO(NFREU, IFIX(10.  /  0X1  +  1.0)) 

NOPT  S  =  M  A  X  0  (NlJPT  S,  20) 

DO  bOOO  1  =  1, NOPTS 
XARR(i)  =  FREU(l) 

Y  APR ( I )  =  PS(I) 

6000  CONTINUE 

NPP I  =  IFIX(1.0  /  UX1) 

NPP 1 -M A  X  0 ( WPP 1,2) 

CALL  OVERLAY (1HK, 3, 1 ) 

7000  CONTINUE 

C******  ***************************** 

9000  RETURN 

9001  FORMATdX,/,  1 0 X , 46H* * N ARN1 NG* *  ERROR  TERM  FROM  SUBROUTINE  P1FREQ  , 

A  abH  ( OEF I NEO  bY  SUBROUTINE)  IS,  15,  //) 

900a  FORMAT  d X , / / / , 20 X , 27 HT EN  PREDOMINANT  FREQUENCIES, 

A  //,5X, 10F10.5,/) 

9003  FUKMAT (312) 

9009  FORMATdX,////,  10X,23HPLOT  OPTION  FOR  PAWS  IS ,  5X  ,  1 5 , 2X ,  '  L  AGOPT  =  • 

*  ,15) 

90  05  FORMAT (1  HI ,///, 30X, 54 (1 H* ) , / , 2 0 X , 1 H* , 53 X , 1 H * , / , 20 X , 1 H* , 1 5X , 

*  23HP0NER  SPECTRUM  ANAL YSI S, 1 4X , 1 H* , / , 20X , 1 H* , 52 X , 

*  1H*,/,20X,54(1H*)) 

END 
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SUBROUT  X NE  PAWOUT 

C - 

C 

C  *  SUBROUTINE  PAYOUT  PKINTS  THE  OUTPUT  FROM  SUBkUUTINE  PAWS 
C 

C - - - FALCON  R&D  (ME  A  b-5-79 

COMMON  /  CONTRL  /  LUIN  ,  LUOUT  ,  LUOATA  ,  LORON  ,  LUREPT 
COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACri 
COMMON  /  CONTRL  /  USRNAM  ,  PUHPOS 
INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,  PURPOS(b) 

COMMON  /  PAwCOM  /  XYMV(b),  ACV(lb02),  FRE0(601),  PS(lb02), 

*  XCOV  ( 1  602)  ,  XSPECT (lb02) ,  AMPHAS  (1602)  , 

*  XFER(lb02),  COHER(eOl) 

LEVEL  2,  LU I N, LUOUT, LUDAT A, LORUN, LUREPT, IFSVRP, ONLINE, IANOPT 

*  ,  0PTI0N,USRNAM,PURPUS,M,LAG,F1N0 

*  ,  XYMV, ACV,FREU, PS, XCOV, XSPECT, AMPHAS, XFER,COHER 

LEVEL  2 , IFB ACH 
COMMON  /PCOMM/  LAG  ,  FIND 
COMMON  /  LAGCOM  /  M 

*  PRINT  HEADING,  MEAN,  VARIANCE 

WRITE  (LUOUT,  9001) 

WRITE  (LUOUT,  9002)  XYMV(l),  XYMV(2) 

C 

C  *  PRINT  LAG  AND  AUTOCOVARIANCE  (LAG)  TABLE 
C  FREUUENC Y  AND  POWER  SPECTRUM  TABLE 

C 

WRITE  (LUOUT,  9003) 

LAG  =  0 

WHITE  (LUOUT,  9004)  LAG,  XYMV(2),  FREQ(l),  PS(1) 

00  1000  LAG=1,M 
FIND  =  LAG  +  1 

WRITE  (LUOUT,  9004)  LAG,  ACV (LAG) ,  FREQ  (FIND) ,  PS(FIND) 

1000  CONTINUE 

WRITE  (LUOUT,  9005) 

9000  RETURN 

9001  FORMAT  (INI,  //,  20X,  34HRESULTS  OF  POWER  SPECTRUM  ANALYSIS,  //  ) 

9002  FORMAT  (IX,  10X,  1 3HWA VEFORM  MEAN,  F15.4,  /,  11X, 

*  1 7HWAVEF  ORM  VARIANCE,  F15.4,  /) 

9003  FORMAT  (IX,  //,  17X,  20 H AUTOCOVARIANCE  TABLE,  10X, 

*  14HP0WER  SPECTRUM,  //,  17X,  3HLAG,  3X,  1 4H AUTOCOV ARI ANCE , 

*  5X,  9HFRE0UENCY,  3X,  14hPOWER  SPECTRUM,  /) 

9004  FORMAT  (10X,  110,  F15.5,  5X,  F10.5,  F15.5) 

9005  FORMAT  (IX,  ///,  25X,  30HEND  OF  POWER  SPECTRUM  ANALYSIS,  ////) 

END 
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SUBROUT  I NE  RUKSRF (A,B,k,N) 


*  SUBROUTINE  ROrtSRT  SORTS  RUNS  OF  A  MATRIX  A  ACCORDING  TO  DESCENDING 
ORDER  UF  MATRIX  B.  STORE  SORTED  A  INTO  R. 


DIMENSION  *  A (81 0) ,  B(810),  R(lb20) 

COMMON  / SKTCOM/  1 , 1 2 ,  I  SDR T , KS A VE , S A VE R , I R  ,  I  A 
LEVEL  2,  A, B,R, I, N, 12, ISORT, RSAVE, SAVER, IR,1A 

*  MOVE  SORTING  KEY  VECTOR  TO  FIRST  COLUMN  OF  OUTPUT  MATRIX  AND 

BUILD  ORIGINAL  SEQUENCE  LIST  IN  SECOND  COLUMN 

DU  10  00  J  =  1  ,  N 
RCI)  =  B(I) 

12  =  I  +  N 
R  ( 1 2  )  =  1 
1000  CONTINUE 

*  SORT  ELEMENTS  IN  SORTING  KEY  VECTOR  (ORIGINAL  SEQUENCE  LIST  IS 

RESEQUENCED  ACCORDINGLY ) 

2000  CONTINUE 
I  SORT  =  0 
DO  0000  I  =  2  ,  N 

IF  ( R  (  I )  -  R  ( 1-1 ) )  4000,  4000,  3000 
3000  CONTINUE 

ISORT  =  ISORT  +  1 
RSAVE  =  R  ( 1 ) 

R ( I )  =  k(I-l) 

R ( I - 1 J  =  RSAVE 
12  =  I  +  N 

SAVER  =  R ( I 2 ) 

R( 12)  =  R ( 12- 1 ) 

R ( 12-1 )  =  SAVER 
4000  CONTINUE 

IF  (ISURT)  2000,  5000,  2000 

*  MOVE  ROaS  FROM  MATRIX  A  TO  MATRIX  R  (NUMBER  IN  SECOND  COLUMN  OF  R 

REPRESENTS  ROw  NUMBER  OF  MATRIX  A  TO  BE  MOVED) 

5000  CONTINUE 

DO  8000.  I  =  1  ,N 

*  MOVE  ELEMENT  TO  OUTPUT  MATRIX 

12  =  I  +  N 

R  ( I  )  =  A  (  R  ( 1 2 )  ) 

8000  CONTINUE 
RETURN 
END 
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OVERLAY (l\,  i,  1  ) 

PROGRAM  SCklbl 

c  — * - - - - 

c 

C  *  SCR1B1  PLOTS  GRAPHS  FOk  MUZZLEDEV ICE  DATA  IN  PArtS 
C  USING  CALCOMP  PLOTTING  PACKAGE 

t 

- - 

COMMON  /  SCklbl  /  XAkRAY (100)  ,  YAkkAY (100)  ,  DX1  ,  NDPTS  ,  NPPI 
REAL  OX  ,  XSCALE  ,  YSCALE  ,XLNGIH. 

REAL  YLNGTH  ,  XOR  ,  YOH  ,  XBASE  ,  YBASE  ,  FACTOR 

INTEGER  I  SENT (10)  ,  L ABEL  C A )  ,  1UNIT 

COMMON  /  LABEL  /  KTB(BA)  ,  TST  ART  ,  I 
LE  VEL  2,  KTb,  TST  ART ,  1 

C 

DATA  LABEL/ 'BOOTS' , 'B390' , 'X3983' , 'M392  DATA'/ 

C 

C  TES1  ARGUMENTS  FOR  ERRORS 

C 

IF  (NDPTS  .GT.  100)  CALL  PERROR  (8P) 

DEFINE  VARIABLES  USED  BY  PLOTTING  KUUTINES 

XLNGTH=1 1 . 

YLNGTH=8.5 
XSIZE=8.0 
YS1ZE=6.0 
XbASE=P.O 
YBASE  =  1 .0 

CALL  FIXSCA (XARRA Y (1 ) , NDPTS, XS1 ZE , XSCALE , XOR , XMAX , DX ) 

CALL  FIXSCA (YARRAY (1 ) , NDPTS, YSIZE, YSCALE, YOR, YMAX,DY) 

XLNGTH=1 1 . 

YLNGTH=8.5 
XS1ZE=8.0 
YSIZE=6.0 
XBASE=2.0 
YBASE  =  1 .0 

CALL  FIXSCA (X  ARRAY (1 ) , NDPT S , XS I ZE , X SC ALE , X OR , XM A  X , D X ) 

CALL  FIXSCA (YARRAY (1) , NDPTS, YSIZE, YSCALE, YOR, YMAX,DY) 

FACTOR  =  1. 

I UN  I T  =  10 

BRL  PLOTTING  ROUTINES  (INITIAL) 

CALL  PL TBEG(XLNGTH, YLNGTH, FACTOR, I  UN IT, LABEL) 

CALL  PLTSCA (XBASE, YBASE, XOR, YOR , XSCALE , YSCALE ) 

PLOT  THE  CURVES 

CALL  PLTDTS(1, 0, XAkRAY, YARRAY, NDPTS, 0) 

PLOT  AXES  AND  LABELS 

MODE  =  A 
ONE  =  1. 

CALL  PLTAXS(DX,DY, XOR, XMAX, YOR, YM AX, MODE) 

CALL  L ABEL A (OX, DY, XOR, XMAX, YOR, YMAX , ONE, ONE) 
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YMAX1=YmAX+1 ,U*YSCALt 
YMAX1=YMAX+1 .0*YSCALL 
ENCODE  (  100,  130,  1  SENT ) 
130  FORMAT  (IX,  2412,  '>') 

CALL  PLTSYM(0. 15, ISENT, 0. 
CALL  PLTSYM(0. lb, ISENT,  0. 

CALL  PLTPGE 

RETURN 

END 
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OVERLAY (E  ,  4 , 0  ) 

PROGRAM  G1RAFE 

C - 

c 

C  *  SUHKOUT 1  ME  GIRAFE  DOES  A  FORWARD  AND  BACKWARD  INTERPOLATION 
C  FROM  EACH  POllvl  IN  THE  STRING  AND  CALCULATES  THE 

C  ERROR  FOR  THE  EXTRAPOLATED  POINT 

C 

C - FALCON  R&U  fa/28/79  NEA 

$  COMMON  /  SSDATA  /  NSSPTS,  SDTIME,  SUBSTK (5000 ) ,  SSTIME(5000) 

$  DIMENSION  ERR(20, 2),  WKAREA  (50)  ,  A(10,4),  8(10,1),  C(10) 

$  LEVEL  2,  NSSPTS,  SDTIME,  SUBSTR,  SSTIME 

$  LEVEL  2,  NL1M, 1ST, LLST, INDEX, IND, XI, IrtK,C,YI,YF,YB, ERR 

$  COMMON  /  CMGIRF  /  NL I M , 1  ST , LLST , I NDE X , I NO , X I , I WK , C , Y 1 , Y F , Y B , ERR 

C 

C  *  READ  NUMBER  OF  DATA  POINTS  FOR  INTERPOLATION  FIT 
C 

$  READ  (5,  9004 )  NBATS 

C 

C  *  CALCULATE  NUMBER  OF  EXTRAPOLATION  POINTS 
C 

$  NLIM  =  NSSPTS  -  2 * N B A T S 

*  NCOEF  =  4 

$  NDIM  =  10 

$  NONE  =  1 

$  NDGT  =  7 

C 

C  *  LOOP  FOR  FORWARD,  BACKWARD  INTERPOLATION 
C  INTERMEDIATE  EXTRAPOLATION 

C 

$  DO  1000  IST=1,NLIM 

C 

C  *  FORWARD  EXTRAPOLATION 
C 

$  LLST  =  1ST 

$  DO  2000  INDEX  =  1, NBATS 

$  IND  =  LLST  +  INDEX  -  1 

$  XI  =  SSTIME(IND) 

$  A ( I NDE X ,  1)  =  X  I  **3 

$  A  ( INDEX  ,  2 )  =  X I  **2 

$  A ( INDEX, 3)  =  XI 

$  A(INDEX,4)  =  1.0 

$  B(INDEX,1)  =  SUBSTR (IND) 


$2000  CONTINUE 
C 

C  *  LEAST  SUUARES  FIT 
C 

$  DO  2500  I WK= 1 ,50 

$  WKAREA (IWK)  =  0.0 

$2500  CONTINUE 

$  CALL  LLSOAR  (A,  B,  NBATS,  NCOEF,  NONE,  NDIM,  NDIM,  NDGT,  WKAREA, 

$  *  IER) 

$  DO  3000  IND=1,4 

$  C(IND)  =  B ( I ND ,  1) 

$3000  CONTINUE 

$  XI  =  SST I  ME ( I  ST  +  NBATS) 

$  YI  =  SUBSTR (1ST  ♦  NBATS) 
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$  YF  =  C(1)*XI**3  +  C(2)*X1**2  +  C(3)*XI  +  C(4) 

4  WRITE  (6,  9001)  (C(I),  1=1, a) 

4  00  3500  INDEX  =  1 ,  NBAT S 

4  INI)  =  LLST  +  INDEX  -  1 

$  XI  =  SSTIME(INU) 

s  ri  =  cu)*xi**3  +  c(2)*xi**2  +  c(3)*xi  +  c(4) 

4  WRITE (6,9006)  XI  ,  SUBSTR(IND)  ,  YI 

$3500  CONTINUE 

$  WRITE  (fa,  9007)  XI,  YI,  YE 

C 

C  *  BACKWARD  EXTRAPOLATION 
C 

$  LLST  =  1ST  +  NBATS  +  1 

$  DU  4000  INDEX  =  1, NBATS 

$  INU  =  LLST  +  INDEX  -  1 

$  XI  =  SSTIME(IND) 

4  AlINDEX, 1)  =  X  I  *  *  3 

4  A(lN D£X»2)  =  X  I  **2 

$  A ( I NDEX , 3)  =  XI 

5  A (INDEX, 4)  =  1.0 

$  B ( INDEX , 1 )  =  SUBSTR (  IND) 

S4000  CONTINUE 
C 

C  *  LEAST  SQUARES  FIT 
C 

$  DO  4500  I WK  =  1 , 5 0 

$  WK  ARE  A ( I WK )  =  0.0 

$4500  CONTINUE 

$  CALL  LLSUAR  (A,  B,  NBATS,  NCOEF ,  NONE,  NDIM,  NDIM,  NOGT  ,  WK  ARE  A , 

$  *  IER) 

4  00  5000  I N D  = 1 ,4 

4  C (IND)  =  b (IND, 1 ) 

$5000  CONTINUE 

$  XI  =  SSTIMECIST  +  NBATS) 

$  YI  =  SUBST  R ( 1ST  +  NBATS) 

$  YB  =  C(1)*XI**3  +  C(2)*XI**2  +  C 1 3 ) *  X I  +  C(4) 

$  WRITE  (fa,  9002)  (C(I),  1=1,4) 

5  WRITE  (fa,  9007)  XI,  YI,  YB 

$  DO  5500  INDEX  =  1, NBATS 

4  IND  =  LLST  +  INDEX  -  1 

$  XI  =  SSTIME(IND) 

5  YI  =  C ( 1 ) *  X 1 **3  +  C(2)*XI**2  +  C ( 3 ) *  X I  ♦  C(4) 

$  WRITE(b,900fa)  XI  ,  SUBS  1 R ( IND)  ,  YI 

45500  CONTINUE 


*  CALCULATE  ERROR 


ERR  ( I  ST ,  1 )  =  SSTIMEUST  NBATS) 

ERR(IST ,2)  =  YF  -  YB 
WRITE  (b,  9003)  ( ERR ( I  ST , I ) ,  1  =  1,2) 

41000  CONTINUE 
9000  RETURN 

$9001  FORMAT  (IX,  //,  21HFORWARD  EXTRAPOLATION,  /, 

$  *  2faHLE AST  SQUARES  COEFF I C 1  ENTS ,  4E14.7) 

49002  FORMAT  (IX,  /,  22HBACKwARD  EXTRAPOLATION,  /  , 

S  *  26HLEAST  SQUARES  COEFFICIENTS,  4E14.7) 

$9003  FORMAT  (IX,  /,  4  FIT  IME,  E14.7,  5X,  22HF0RW  ARD-B  ACKWARD 


ERROR, 
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4  *  fc 1 4. 7 , ///) 

4^004  FUKMAT  (12) 

490  0b  FUKMAT  (IX,  37HLST  SDK  DATA:  INDEX,  X3  ,  X2»  X»  1,  Y  , lb»5E14.  ) 

5,900b  FUkMAT  ( 1  X ,  4 H T  I ME ,  2 X  ,  1 1  4 . 7 , 2 X  ,  bH V  ALUE ,  2X  ,  fc  1  4 . 7 , 2 X  , 
j  *  tiHESTlMATE,2X,E14.7) 

49007  FUkMAT  (IX,/,  4HT IME , 5X , E 1 4 . 7 , 1  OX , 5H V ALUE , bX , E 1 4 . 7 , 1  OX , 

4  *  ttHEbTIMATE,5X,fcl4.7,/) 

END 
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OVERLAY (F ,5,0) 
PROGRAM  thefft 


C 


*  SUBROUTINE  THtFFT  ANALYSIS  IS  FOk  THE  FIRST  OR  THIRD  REGIME 
POLY NUMl AL  FIT,  FAST  FOUklfck  TRANSFORM 
REMAINDER  AFTER  SUBTRACT ION  OF  FITS  IS  NOISE 


CONTRL  /  LUIN  ,  LUOOT  ,  LODATA  ,  LURUN  ,  LUREPT 
CONTRL  /  1FSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH 
CONTRL  /  OSRNAM  ,  PURPOS 

SSDATA  /  NSSP1S,  SOT  I  ME ,  SOBSTR (5000) ,  SSTIME(5000) 

SOATA  /  NSTPTS  ,  DTIMt  ,  STRING15000)  ,  STIMEC5000) 
REPOAT  /  REPSTR (5000)  ,  IFREP  ,  NDREP 
FREUS  /  FREQEN(20),  NFR 
CMLOOT  /  AC(20),  BS (20 ) ,  INDEX ( 1 50Q ) , 

CSTR11500),  OST  k ( 1500) ,  PER  ( 1500) , 

RN01SE (1500) ,  1,  AO,  OMEGA,  PI, 

J,  K,  SUMER,  AVEER,  SIGMA,  PL0T0P(2), 

NBIN,  CELLS ( 1 0 ) ,  COMP(IO),  CS,  IDF,  0, 

1 1  NO , SUMD ,  AVEO,  SIGD,  SIGNDD ,  STAR, 

TIME,  1  REP  ,  1ST  ART ,  ISTOP,  11,  12 
SCRIB  /  X ARR  (5000)  ,  YAkRl(5000)  ,  YARR2(5000), 

Y  ARR3 (5000)  ,  NARR  ,  DELTAX  ,  NP01NT  ,  NPPI 

POLYCM  /  BS1 R ( 1 50 0 ) ,  MAXDEG,  POLY(1500),  COEF(IO),  NCOEF 
DER1VS  /  STRDOT (5000) ,  STRDD(5000),  PLT0P1 

NSSPTS,  SDTIME,  SOBSTR,  SSTIME,  FREOEN,  NFR, 

AC,  HS,  MAXDEG,  COEF,  BSTR,  CSTR,  DSTR,  PER, 

RNOISE,  POLY,  I,  AO,  OMEGA,  PI,  J,  K,  SUMER, 

AVEER,  SIGMA,  NCOEF,  PLOTOP,  INDEX, 

NBIN,  CELLS,  COMP,  CS,  IDF,  Q,  NDREP, 

NST PTS, DTI  ME, STRING, ST  1  ME, REPSTR, IFREP, 

STRDOT,  STROD,  PL TOP  1 ,  1 1 NO,  SUMD,  AVED,  SIGD, 

SIGNDD,  STAR, T IME, IkEP, 1ST  ART, ISTOP, 11 , 12 
LEVEL  2, LUIN, LUOUT, LUD AT A, LURUN, LUREPT, I FS VHP, ONLINE 
LEVEL  2, IANOPT, OPTION, I FB ACh , USkN AM , PURPOS 
C 

INTEGER  PLOTOP  ,  PLTOP1 
EXTERNAL  NORM 
C 

c*********************************** 

c 

C  *  SET  UP  CONSTANTS 
C 

NOMA  X  =  1500 
NOPTS  =  NSSPTS 
PI  =  3.1015926 
NCOEF  =  0 
NBIN  =  6 
IDF  =  0 
MAXDEG  =  0 

IF  (NSSPTS  .GT.  NDMAX)  CALL  PERROR  (81) 

C 

C  *  WRITE  PAGE  HEADING  FOR  THEFFT 
C 
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WRITE (6, 90  IB) 


C 

C  *  PLUT0P(2)  =  PLOT  OPTIONS  F  Ok  DATA  SUBSET  ANALYSIS 

=  1  ON 
=  0  OFF 

(1)  s  3  PLOTS  OF  DATA,  POLYNOMIAL  FIT,  PERIUD1C  FIT 
AND  NOISE 

(2)  =  1  PLOT  OF  DATA  SET  AF 1  EH  NOISE  HAS  BEEN  REMOVED 

HEAD (5, 9002)  PLOTOP 
WR I TE (fc, 9009)  PLOTOP 
C 

c 

c 

C  *  BSTH  IS  ORIGINAL  STRING 
C 

DO  1000  1  =  1 , NOPTS 
BSTRll)  -  SUBSTR  ( I ) 

1000  CONTINUE 
C 

c***************************^^^^^^^^ 

c 

C  *  FOURTH  DEGREE  ORTHOGONAL  POLYNOMIAL  FIT  UN  BSTR1  IN  POLY 
C  CSTR  IS  ORIGINAL  MINUS  POLYNOMIAL,  WHICH  IS  THE  LEVELED 
C  CURVE  OF  THE  ORIGINAL. 

C 

CALL  OVERLA Y ( 1 HH, 5, 1 ) 

WKITElb, 901-1)  (COEF(I)  ,  I  s  l,NCOEF) 

DO  2000  I  =  1 , NOPTS 

CSTR  (  I  )  =  BSTR(I)  -  'POLY  Cl) 

SUBSTR(l)  =  CSTR(I) 

2000  CONTINUE 

c*****************************aa^^aa 

c 

c  *  OBTAIN  10  MAJOR  FREUENC1ES  OF  CSTR  IN  FREQEN  BY 
C  SR  PAWS  -  POWER  SPECTRUM. 

C 

CALL  OVERLAY (1HI,5,2) 

C 

C  *  CALCULATE  COEFFICIENTS  OF  COSINES,  SINES  OF  MAJOR  FREQUENCIES 
C 

AO  =  0. 

C 

C  *  CONSTANT  AO 
C 

DO  2750  1  =  1 , NFR 

IF  (FREUEN(l)  .NE.  0)  GO  TO  2750 
AO  =  0. 

GO  TO  3500 
2750  CONTINUE 

DO  3000  J  =  1 , NDPTS 

AO  =  AO  +  2.  *  CSTR(J) 

3000  CUNTINUE 
3500  CONTINUE 
C 

C  *  COEFFICIENTS  AC(K)  AND  BSIK)  FOR  COSINE  AND  SINE  OF 
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C  FREQUENCY  FREQEN(K) 

C 

WkITE(b,90l9) 

WRITE(b,90l3) 

UU  4000  I  =  1,NFR 
AC ( 1 )  =  0. 

BS  Cl)  =  0. 

DO  4500  J  =  1  , N DPT S 

OMEGA  =  2.  *  PI  *  (FLOAT ( J- 1 )  /  FLOAT (NDPTS) )  *  (FREQEN(l) 

A  *  SDTIME  *  NDPTS) 

AC ( 1 )  =  AC ( I )  +  CSTR(J)  *  COS ( OMEGA) 

bS  ( I )  =  bS  ( I )  +  CSTR(J)  *  SIN  IOMEGA) 

4500  CONTINUE 

AC ( I )  =  2.  *  A  C ( I ) 
bS ( I )  =  2.  *  BS(I) 

WR I TE (b , 90 1 4 )  I,  FREQEN(I),  AC(I),BS(I) 

4000  CONTINUE 

C********************************* 

c 

C  *  RECONSTRUCT  PERIODIC  WAVEFORM  FROM  MAJOR 
C  FREQUENCIES. 

C  CALCULATE  PER(I)  FOR  POINT  I. 

C 

DO  5000  J  =  1, NDPTS 
PER  ( J )  =  0. 

DO  7000  K  =  1,NFR 

OMEGA  =  2.  *  PI  *  (FLOAT (J-l )  /  FLOAT (NDPTS) )  *  (FREQEN(K) 

A  *  SDTIME  *  NDPTS) 

PERU)  =  PER  ( J )  +  A  C  ( K )  *  COS(OMEGA)  +  BS(K)  *  SIN  (OMEGA) 

7000  CONTINUE 

PER  ( J  )  =  PER  ( J  )  /  FLOAT (NDPT S) 

PER  ( J )  =  PER  ( J )  +  AO  /  (2.  *  FLOAT (NDPTS) ) 

5000  CONTINUE 

C*********************************** 

c 

C  *  NOISE  ARRAY  IS  LEVELED  CURVE  MINUS  PERIODIC  FIT  TO  IT. 

C 

DO  bOOO  I  =  1, NDPTS 

DSTR(l)  =  CSTR(I)  -  PER  ( I ) 

SUBSTR ( I )  =  DSTR(I) 

SUMER  =  SUMER  +  OSTR(I) 
bOOO  CONTINUE 

AVEER  =  SUMER  /  NDPTS 
C 

C*********************************** 

c 

C  *  PRINT  ARRAYS  OF  ORIGINAL  DATA,  ORIGINAL  ESTIMATE,  DIFFERENCE, 

C  POLYNOMIAL  FIT,  DIFFERENCE,  PERIODIC  FIT,  AND  DIFFERENCE  (NOISE) 

C 

WRITE(b,901b) 

WRITE (b, 9003) 

DO  0000  I  =  1, NDPTS 

WRITE(b,9004)  I,  SSTIME(I),  BSTR(I),  POLY(I), 

*  CSTR(I),  PER ( I ) ,  DSTR(I) 

8000  CONTINUE 

c*********************************** 

C 
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C  *  CALCULATE  S1ANUAR0  DEVIATION  UF  NOISE 
C 

SIGMA  =  U. 

DU  0300  1  =  1, HURTS 

SIGMA  =  SIGMA  +  ( DST  R  (  I )  -  AVEER)  **3. 

8300  CONTINUE 

SIGMA  =  SORT (SIGMA  /  FLOAT (NDRTS) ) 

*********************************** 

*  NORMALIZE  FIRST  DERIVATIVE 

IINO  =  (SSTIME(l)  -  STIME(l))  /  DT 1  ME  +  1 
IF  (STIME(IINU)  .NE.  SSTIME(l))  CALL  PERR0R(6O) 

*  AVERAGE  IN  AVED,  STANDARD  DEVIATION  IN  SIGD 

SUMD  =  0.0 

DO  7000  I  =  1 , NSSPT S 
J  =  1IND  +  I  -  1 
SUMD  =  STRDOT (J)  +  SUMD 
7000  CONTINUE 

AVED  =  SUMD  /  FLOAT (NSSPTS) 

SIGD  =  0. 

DO  7600  I  =  1, NSSPTS 
J  =  IINO  +1-1 

SIGD  =  SIGD  +  (ST ROOT ( J )  -  AVED)  **  3 
7600  CONTINUE 

SIGD  =  SORT (SIGD  /  FLOAT (NSSPTS)  ) 

*************  ********************** 

*  NORMALIZE  NOISE  AND  ANALYZE  IT 

WRITE (6, 9015) 

DO  8050  I  =  1 »  NDP1 S 

RNOISE(I)  =  (USTR(l)  -  AVEER)  /  SIGMA 
8050  CONTINUE 

CALL  GFIT (NORM, NB IN, RNOISE, NSSPTS, CELLS, COMP, CS, IDF, U, IER) 

IF  (IER  .GT.  0)  CALL  PERR0R(85) 

WRITE(6,9005)  AVEER  ,  SIGMA 

WRITE (6, 9007)  (CELLS(I),I  =  1,NBIN) 

WRITE (  6 ,  9008)  (COMP(I),  I  =  1,NBIN) 

WRITE (6, 9006)  IDF  ,  CS  ,  0 
WRITE  (6,9017) 

J  =  0 

DO  8000  I  =  1, NDRTS 
STAR  =  1H 

IF  (RNOISE (I )  .GT.  3.0)  STAR  =  1H* 

SIGNUD  =  ( STRDOT ( 1 )  -  AVED)  /  SIGD 

IF  (RNOISE (I)  .LT.  3.5  .AND.  SIGNDD  .LT.  3.0)  GO  TO  8000 
WRITE (b, 90 18)  I,SSTIME(I) ,DSTR (1) , RNOISE (I) , STAR, SIGNDD 
8000  CONTINUE 


*  CALCULATE  INDEX  IN  STRING  CORRESPONDING  TO  START  OF  ARRAY  SU6STR 
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TIME.  =  SSTIMEU) 

IF  (SDTIME  .NE.  UTlMfc)  CALL  PERROR  (63) 

I1ND  =  (TIME  -  ST IME ( 1 ) )  /  DTImE  +  I 

IF  (ST  IME  1 1  INI)  J  .NE.  SSTIMt(l))  CALL  PERROR  (64) 

*  REPLACE  SUBSTK  BY  ITS  (HEORETICAL  APPRUX 1MAT ION  UF  A 

POLYNOMIAL  PLUS  A  PERIODIC. 

DU  5500  I  =  1,NSSPTS 
IREP  =  1IM0  +1-1 
REPSTR ( I  REP )  =  POLY(I)  +  PER(I) 

500  CONTINUE 

*********************************** 

*  SET  IF  REP  TU  INDICATE  PROGRESS. 

IF  (IF  REP  .ED.  2)  IFREP  =  3 

IF  (IFREP  .  E  !J .  1)  IFREP  =  2 

IF  (IFREP  .EO.  0)  IFREP  =  1 


IF  ( PLOT  OP ( 1 )  .EU.  0)  GO  TO  6900 

PLOT  DATA  SUBSET ,  POLYNOMIAL  FIT 

DO  8500  1  =  1  , NUP T S 
YARRl(I)  =  BSTR(l) 

YARR2(I)  =  POLY(I) 

YARR3(I)  =  0 
XARR(I)  =  SSTIME(l) 

8500  CONTINUE 

DELTAX  =  SDTIME 
NPOINT  =  NDPTS 
NARK  =  2 
NPP I  =  40 

CALL  0VERLAY(1HJ,5,3) 

PLOT  DATA  SUBSET  -  POLYNOMIAL  FIT  ,  PERIODIC  FIT 

NARR  =  2 

DO  8750  I  =  1, NDPTS 
YARR1 (I)  =  CSTR(l) 

YARR2(I)  =  PER ( I ) 

YARR3 ( I )  =  0. 

8750  CONTINUE 

CALL  OVERLAY (1HJ, 5,3) 

PLOT  NOISE  =  DATA  SUBSET  -  POLYNOMIAL  FIT  - 
PERIODIC  FIT 

NARR  =  1 

DO  6800  I  =  1, NDPTS 
YARklll)  =  DSTR(I) 

Y  AKR2 ( I )  =  0. 

Y  ARR3 ( I )  =  0. 

PROGRAM  LISTING  (Continued) 


75 


oooo  oono  ooo 


8800  CONTINUE 

CALL  OVERLAY (1HJ, 5, 31 

C********************************* 

89  0  u  CONTINUE 

IF  ( IFRER  .NE.  NUhEP)  GU  10  9000 
IF  ( PLOT  OP (P )  . EQ .  0)  GO  TO  8990 

*  PLOT  DATA  RECORD  WITH  NOISE  REMOVED 

IF  (IANOPT  .EO.  PI  GO  TO  89P0 

THEFFT  CALLED  FROM  KITTY8 
PLOT  ONLY  CURRENT  DATA  SET 

ISTART  =  I1NU 

ISTOP  =  IIND  +  NSSPTS  -  1 
NPUINT  =  NSSPTS 
GO  TO  8950 
89P0  CONTINUE 


THEFFT  CALLED  FROM  COMPLT 
PLOT  ENTIRE  DATA  SET 


ISTART  =  1 
ISTOP  =  NSTPTS 
NP01NT  =  NSTPTS 
8950  CONTINUE 
NARR  =  1 

DO  6500  I  =  ISTART,  ISTOP 
XARR(I)  =  STIME(l) 
YARR1 (I)  =  REPSTR(I) 
YARRPC1)  =  0. 

YARR3CI)  =  0. 

6500  CONTINUE 

CALL  0VERLAY(1HJ,5,3) 
8990  CONTINUE 


C 

9000  RETURN 
900P  FORMAT (PIP) 

9003  FORMAT  ( 1  X , / / / ,  1 X , 5H 1NDE X ,  5X , 4HT IME , 1 X , 5X , 6HSTR 1 NG , 5X , 

*  4X , 8HPOL Y  FIT,4X,2X, 13HLE VELEU  CURVE, IX, 

*  4X,8HPERI0DIC,4X,5X,5HN01SE,//) 

9004  FORMAT  ( 1  X ,  1 5 ,  F 1  0 . 5 ,  7(2X,E14.7)) 

9005  FURMAT (IX, ///, 10X, 1 3HMEAN  OF  NO I SE , 1 7 X , F 1 0 . 3 , / , 

A  10X, 27 H STANDARD  DEVIATION  OF  NO  I  SE , 3X , F 1 0 . 3 , / ) 

9006  FORMAT  ( 1  HO ,  1  0 X ,  1 8HDEGREES  OF  F REEDOM , 3 1  X , 1 1 0 , / , 

*  15X, P1HCH1-SOUARED  S T A T 1 S T 1C , 28X , F 1 0 . 3 , / , 

*  15X, 21HCH1-S0UARED  STATISTIC,/, 

*  15X,  5X, 39HEXCEEDING  CS  IF  NULL  HYPOTHESIS  IS  TRUE , 5X , F 1 0 . 3 

*  //) 

9007  FORMAT  (10X,8HCELLS  =  ,10F10.3) 

9008  FORMAT  (10X,8HC0MP  =  ,10F10.3) 

9009  FORMAT (IX,////, 10X,25HPLOT  OPTIUN  FOR  THEFFT  1S,5X,PI5) 

9011  FORMAT (IX,///, 1 5X , 54HC0EFF IC 1 EN T S  OF  FITTED  POLYNOMIAL  (Cl,  C2,... 
*,  CM,  CO) ,/, 10X,8E15.7, //) 
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901  a  FORMAT (1H1,///,20X,54(1H*),/,20X,1H*,52X,1H*,/,20X,1H*,13X, 

*  27HREGIME  FITTING  AND  AN AL YSIS , 1 2X , 1 H* , / , 20X , 1 H* , 

*  S2X, 1H*,/,20X,54 (1H*) ) 

90  13  FORMAT  (  1  X  ,/////,  32  X  ,  9HF  RE  NUENC  Y  ,  1  bX  ,  bHCOS  I NE  ,  1  bX  ,  4HSINE , 

*  /,3bX,bH(kHZ) , 16X, 1 1HCOEFF 1C  I  ENT , 14X, 

*  1 1HCOEFFICIENT, //) 

9014  FORMAT ( 10X,  1 10,  i ( 1  OX , E 1  4 . 7  )  ) 

901b  FORMAT (1H1,///,20X,54(1H*J,./,20X,1H*,52X,1H*,/,20X, 

*  lrt*, 18X, 17HANALYS1S  OF  NOISE, 17X, 1H*,/,20X, 

*  1H*,52X,1H*,/,20X,S4(1H*)) 

901b  FORMAT (lHl,///,20X,54(lH*),/,20X,lH*,b2X,lH*,/,20X,lH*, 

*  20X, 13HFITTEO  CUk VES , 1 9x , 1 H* , / , 20 X , 1 H* , S2X , 1 h* , 

*  / , 20  X , b4 ( 1 H* ) ) 

901/  FOkMAT (IX, /////, 31X,29HN01SE  GREATER  ThAN  2.b  SIGMAS,/, 

*  faX,bHINDEX,SX,4Hl 1ME, 1 2 X , 5HN0 I SE , b X , 10HNO.  SIGMAS, 

*  5X ,  1 OHGT  3  SIGMA, bX,22HSIG MAS  OF  FIRST  UERIV.,//) 

90  18  FORMAT (6X, Ib,bX,F10.4,bX,F10.4,bX,F10.4, 1 0  X , A  1 , 10X,F10.4) 
9019  FURMAT(1H1,///,20X,54(1HcJ,/,20X,1H*,S2X,1H*,/,20X,1H*,19X, 

*  14HRERI00IC  CURVE, 19X, 1H*,/,20X, 1H*,52X, 

*  1m*,/,20X,54  (1H*)  J 
END 
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non  ooo  o  o  r)  o  ooooooo 


S  U  B  R  0  U  T I  N  E  F  UK  IT  (Nr  M,  ltk) 


*  SUBkOUT 1  NE  F  UR  IT  DOES  A  FOURIER  ANALYSIS  OF  A  TABULATED 

FUNCTION  ANO  KElUKNS  THE  FOUKIEK  COEFFICIENTS. 

*  IBM  SUBROUTINE  PACKAGE 


COMMON  /  LtVELA  /  A(20),  B(20),  FNT(1500) 

COMMON  /FOKCOM/  AN , COEF , CONST , S 1 , C 1 , C , S , J , FNT Z , U 1 , U2 , I , UO , Q 
LEVEL  2,  N,M»IEK,G,U0,U1,U2, I , F NT Z , J , S , C , C 1 , S 1 , CONST , COEF , AN 

LEVEL  2,  FNT,  A,  b 

*  CHECK  FOK  PARAMETER  ERRORS 

IER  =  0 

20  IF  (M)  30,40,40 
30  IER  =  2 
RETURN 

40  IF  (M-N) . 60,60,50 
50  IER  =  1 
RETURN 

*  COMPUTE  AND  PRESET  CONSTANTS 
60  AN  =  N 

COEF  =  2.0  /  (2.0  *  AN  +  1.0) 

CONST  =  3.141593  *  COEF 
SI  =  S I N ( CONST  ) 

Cl  =  COS(CONST) 

C  =  1.0 
S  =  0.0 
J  =  1 

FNTZ  =  FNT(l) 

70  02  =  0.0 
U1  =  0.0 
1  =  2  *  N  +  1 

*  FORM  FOURIER  COEFFICIENTS  RECURSIVELY 

75  UO  =  FNT ( I )  +  2.0  *  C  *  U1  -  U2 
U2  =  U1 
U1  =  UO 
1  =  1-1 
IF  (1-1)  80,80,75 

80  A ( J )  =  COEF  *  (FNTZ  +  C*U1  -  U2) 

8 ( J )  =  COEF  *  S  *  U1 
IF  (J-IM+l))  90,100,100 
90  Q  =  Cl  *  C  -  SI  *  S 
S=C1*S+S1*C 
C  =  0 
J  =  J  +  1 
GO  TO  70 

100  A ( 1 )  =  A ( 1 )  *  0.5 
RETURN 
END 
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o  o  n  d  o 


SUUKOU 1 1 NE  NORM  (X,P) 


*  NORM  CALLS  THt  ImSL  SOBKOUTINt  NH1CH  CALCULATES  Tht 
GAUSSIAN  l)IS  I  KIBU1  ION  FOR  A  VALUE  X 


CALL  MDNOK  U,P) 

RETURN 

END 
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n  o  o  o  o  o 


OVERLAY (H, 5,  1  ) 
PRObRAM  PULYOR 


*  REGRESSION  ANALYSIS  ON  Y  AlIH  RESPECT  TO  X, 

POLYNOMIAL  (DEGREE  NCUEF-1),  COEFFICIENTS  IN 
COEF  f  PREDICTION  IN  YPRED. 


C - - - 

COMMON  /  SSDATA  /  NSSPTS  ,  SDTIME  ,  SUBSTR(5000)  ,  SSTIME(5000) 
COMMON  /  POLYCM  /  Y(1500)  ,  MAXDEG  ,  YPRE0(1500)  ,  COEF(IO)  , NCOEF 
COMMON  /  CMPOLY  /  I , X Y fc ( 1 50 0 , 7 ) , MDP ( 3 ) , NDMA X , NPT , RSQ , AL BP ( 2 ) , 

*  ANOVA( 13) fB(7 , 12) , IB#PRED(1500»6) ,  X (  1 500 ) , I P , WK 
DOUBLE  PRECISION  AK(bOOO) 

LEVEL  2,  Y,  YPRED,  MAXDEG,  COEF,  NCOEF ,  NSSPTS,  SuTIME,  SUBSTR, 

*  SSTIME  ,  XYa,  MDP,  NDMA  X ,  NPT,  RSQ,  ALBP ,  ANOVA,  B,  IP, 

*  PREU,  X,  I,  WK 
LEVEL  2, IB 

C 

DO  1000  I  =  1, NSSPTS 
X(I)  =  SSTIME(I) 

X  Y  W ( 1 ,  1)  =  X  ( I ) 

XYW(I,2)  =  Y(l) 

X Yw ( I , 3)  =  1. 

1000  CONTINUE 

NDMAX  =  1500 
RSQ  =  100. 

IB  =  7 
IP  =  1500 
MUP(l)  =  MAXDEG 
MuP ( 3 )  =  1 
ALBP ( 1 )  =  0.05 
ALBP (2 )  =  0.05 
NPT  =  NSSPTS 
C 

C  *  OBTAIN  FITTED  POLYNOMIAL  USING  ORTHOGONAL  POLYNOMIALS 
C 

CALL  RLFOR  (  X  Yin ,  NDMA X ,  NPT  ,  RSQ,  MDP ,  ALBP ,  ANOVA ,  B ,  I B ,  PRED , 

*  I P , aK , I ER) 

IF  ( I  ER  .Ell.  33)  CALL  PERROR  (87) 

IF  ( IER  .GT.  0)  CALL  PERROR  (88) 

NCOEF  =  MDP (2)  +  1 
DO  2000  1  =  1, NCOEF 
COEF (1)  =  B ( I , 2) 

2000  CONTINUE 

DO  3000  I  =  1, NSSPTS 

YPRED (I)  =  COEF (1)*X(1)  +  COEF ( 2 ) *  X ( I ) *  * 2  +  COEF ( 3 ) *  X ( I ) * *3  + 

*  C0EF(4)*X(I)**0  +  COEF (5) 

3000  CONTINUE 

RETURN 

END 
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ooooooonn  no  oonooo 


OVERLAY l 1,5,2) 
program  paws 

C - - - --- 

c 

*  SUBROUTINES  PAWS  CALCULATES  1  HE  POwER  SPECTRUM  FUR  ARRAY  SUBSTR 
USING  SYSTEM  SUBROUTINE  F1FREQ 
CALLED  BY  SUBROUTINE  THEFFT 
PAWS  =  POWER  ANALYSIS  OF  WAVEFORM  SPECTRA 

- FALCON  R&D  N£ A  b-5-7 

COMMON  /  CONTRL  /  LUIN  ,  LUOUT  ,  LUDATA  ,  LURUN  ,  LUREPT, 

*  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH, 

*  USRNAM  ,  PURPOS 
COMMON  /  SSDATA  /  NSSPTS  ,  SDTIME  ,  SUBSTR (5000)  ,  SST IME (5000) 
COMMON  /  PAwCOM  /  XYMV(b),  ACV(1602),  FKEG(801),  PS(lb02), 

*  XC0V(lb02),  XSPECT  ( 1602) ,  AMPHAS ( 1 602) , 

*  XFER ( 1 6o2) »  COHER (801) 

COMMON  /FREQPS/  A(810)  ,  B(810)  ,  R(1620)  ,  NROW 

Common  /freus/  kreqen(20)  ,  nfr 

COMMON  /  SCRIB2  /  XARR(IOO)  ,  YARR(IOO)  ,  DX1  ,  NDPTS  ,  NPP1 
COMMON  /  LAGC.OM  /  M 

COMMON  /  PAwVAK  /  X(1500)  ,  I  NO  C  6 )  ,  XIND(2)  ,  PLOTOP , I , NFREU 

INTEGER  ONLINE  ,  OPTION  ,  USRNAM(2)  ,  PURPOS ( 6 )  ,  PLOTOP 

LEVEL  2,  NSSPTS,  SDTIME,  SUBSTR,  SSTIME,  A,  B,  R,  FREQEN,  NFR 

*  ,  LUIN, LUOUT , LUDATA, LURUN, LUREPT, IFSVRP, ONLINE,  IANOPT 

*  ,  OPT  ION, USRNAM, PURPOS, M, NROW 

*  ,  1FBACH,  XYMV,  ACV,  FREfJ,  PS,  XCOV,XSPECT,  AMPHAS,  XFER,  COHER 

*  ,  X, IND, X1ND, PLOTOP, I, NFREU 

*  WRITE  PAGE  HEADING  FOR  PAwS 
WRITE (6, 9005) 


*  PLOTOP  =  PLOTTING  OPTION  FOR  POWER  SPECTRUM 
=  0  OFF 
=  1  UN 

READ (5, 9003)  PLOTOP 
WRITE (6, 9004)  PLUTOP 
C 

C  *  PARAMETERS 

C  M  =  NUMBER  OF  LAG  TIMES 

C 

M  =  2.  *  (1.  /  SOTIME) 

IF  (M  .GT.  NSSPTS)  M  =  M  /  2 
IF  (M  .GT.  NSSPTS)  M  =  M  /  2 
C 

C  *  INITIALIZE  FTFREU  ARRAYS 
C 

DO  1000  1=1, NSSPTS 
X ( I )  =  SUBSTR  ( I ) 

1000  CONTINUE 

1ND(1)  =  0 

IND (2)  =  NSSPTS 
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I N  0  ( 3  )  =  0 

INI)  (4)  =  M 

1ND(5)  =  o 

INI)  (6)  =  0 

XlNU(l)  =  SuTIME 
XlN0(2)  =  0 

*  CHECK  INPUT  TO  SUBROUTINE  FTFREQ 

IF  (IN0(2)  .LE.  3)  CALL  PERROR(90) 

IF  (1ND(4)  .LE.  2)  CALL  P£RR0R(91) 

IF  (IND14)  .GE.  I N 0 ( 2 ) )  CALL  PERR0R(92) 

*  DO  SPECTRUM  ANALYSIS 

CALL  FTFREG  (X,  I ND ,  XINO,  XYMV,  ACV,  FREQ,  PS, 

*  XCUV,  XSPECT,  AMPHAS,  XFER,  COHER,  1 E  R ) 

IF  (IER  .NE.  0)  WRITE  (LIJUUT,  9001)  IER 

*  PRINT  OUTPUT  FROM  POWER  SPECTRUM 

CALL  PAWOUT 

*********************************** 

*  PUT  TOP  10  MAJOR  FREQUENCIES  (.FREQUENCIES  WITH  LARGEST 

POWER  SPECTRUM)  INTO  FREQEN. 

NFR  =  10 
NKOW  =  M  +  1 

DO  3000  I  =  1 , N R 0 W 
A (I)  =  FREQ  (I) 

B ( I )  =  PS(IJ 
3000  CONTINUE 

*  SORT  FREQUENCIES  INTU  R  BY  DECREASING  POWER 

CALL  ROWSRT (A, B, R, NROW) 

*  PUT  10  LARGEST  (FIRST  10)  INTO  FREQEN 

DO  4000  I  =  1 , NFR 
FREQEN ( 1 )  =  R ( I ) 

4000  CONTINUE 

*  WRITE  10  PREDOMINANT  FREQUENCIES 

wRITE(b,9002)  ( F REUEN ( I ) , I = 1 , NFR ) 

*********************************** 

IF  (PL0TUP  .EU.  0)  GO  TO  9000 

*  PLOT  LOG(POWER)VS  FREQUENCY  (FIRST  25) 

NFREQ  =  M  +  1 
DO  5000  I  =  1, NFREQ 

IF  (PS(I)  .LT.  1)  PS ( I )  =  1. 

PS(1)  =  ALUG 1 0 ( PS ( 1 ) ) 

PROGRAM  LISTING  (Continued) 


82 


5000  CONTINUE 

0X1  =  FREQ(2)  "  FREG(l) 

N  OPTS  =  MI  NO  (NFREU,  IFIXllO,  /  0X1  +  1.)) 

00  6000  I  =  1 , NOPT S 
X  ARR  111  =  FkEU(l) 

YARk(l)  =  PS(I) 

6000  CONTINUE 

NPPI  =  IF IX ( 1  .  /  0X1) 

CALL  SCRIB2 

C*******' ******  ********************** 

9000  RETURN 

9001  FORMAT ( 1 X ,/, 1  OX » 46H** WARNING**  ERROR  TERM  FROM  SUBROUTINE  FTFREU 

A  2bH ( DEF I NEO  BY  SUBROUTINE)  IS,  15,  //) 

9002  FORMAT  ( 1 X , / / / , 20 X , 27 HI EN  PkEOOMINANT  FREGUENCIES, 

A  // »5X, 10F10.5, /) 

9003  FORMAT (12) 

9009  FORMAT (IX,////, 10X, ,23HPL0T  OPTION  FOR  PAWS  IS,5X,I5) 

9005  FORMAT (1H1,///,20X,59(1H*),/,20X,1H*,52X,1H*,/,20X,1H*,15X, 

*  23HP0WER  SPECTRUM  AN AL YS1 S, 1 4X , 1 H* , / , 20 X , 1 H* , 52X , 1 H* , / , 

*  20X, 1H*,52X, 1H*, /,20X,59 (1H*) ) 

END 
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SUHkUU  r 1 Nt  PAwOUT 


c - - - 

c 

C  *  SUBROUTINE  PAwOUT  PRINTS  THE  OUTPUT  FROM  SUBROUTINE  PAwS 
C 

c - - - FALCON  R&D  NEA  6-5-79 

COMMON  /  CONTRL  /  LUIN  ,  LUOUT  ,  LUDATA  ,  LURUN  ,  LUREPT 
COMMON  /  CONTRL  /  IFSVRP  ,  ONLINE  ,  IANOPT  ,  OPTION  ,  IFBACH 
COMMON  /  CONTRL  /  USRNAM  ,  PURPOS 
INTEGER  ONLINE  ,  OPTION  ,  USRNAM (2)  ,  PURPOS(b) 

COMMON  /  PAWCOM  /  XYMV(6),  ACV(lb02),  FREQ(801),  PS(1602), 

*  XCUV  ( 1602) ,  XSPECT  (lb02) ,  AMPHAS ( 1 602) , 

*  XFER(1602),  COHER(BOl) 

LEVEL  2,  LUIN, LUOUT ,LUDAT  A, LURUN, LUREPT, IFSVRP, ONLINE, IANOPT 

*  ,  OPTION, USRNAM, PURPOS, M, LAG, FIND 

*  ,  XYMV, ACV, FREQ, PS, XCOV, XSPECT, AMPHAS, XFER, COHER 
LEVEL  2, IFBACH 

COMMON  / P C 0 M M /  LAG  ,  FIND 
COMMON  /  LAGCOM  /  M 
C 

C  *  PRINT  HEADING,  MEAN,  VARIANCE 
C 

WRITE  (LUOUT,  9001) 

WRITE  (LUOUT,  9002)  XYMV(l),  XYMV(2) 

C 

C  *  PRINT  LAG  AND  AUTOCOVARIANCE  (LAG)  TABLE 
C  FREQUENCY  AND  POWER  SPECTRUM  TABLE 

C 

WRITE  (LUOUT,  9003) 

LAG  =  0 

WRITE  (LUOUT,  9009)  LAG,  XYMV(2),  FREQ(l),  PS(1) 

DO  1000  LAU=1,M 
FIND  =  LAG  +  1 

WRITE  (LUOUT,  9009)  LAG,  ACV (LAG) ,  FREQ (FIND) ,  PS(FIND) 

1000  CONTINUE 

WRITE  (LUOUT,  9005) 

9000  RETURN 

9001  FORMAT  ( 1 H 1  ,  //,  20X,  39HRESULTS  OF  POWER  SPECTRUM  ANALYSIS,  //  ) 

9002  FORMAT  (IX,  10X,  13H WAVEFORM  MEAN,  F15.9,  /,  11X, 

*  17H WAVEFORM  VARIANCE,  F15.9,  /) 

9003  FORMAT  (IX,  //,  17X,  20HAUT0COVAR1 ANCE  TABLE,  10X, 

*  19HP0WER  SPECTRUM,  //,  17X,  3HLAG,  3X,  1 9H AUTOCO V AR I ANCE , 

*  5X ,  9HFREQUENCY,  3X,  19HP0WER  SPECTRUM,  /) 

9009  FORMAT  (10X,  110,  F15.5,  5X,  F10.5,  F15.5) 

9005  FORMAT  (IX,  ///,  25X,  30HEND  OF  POWER  SPECTRUM  ANALYSIS,  ////) 

END 
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o  o  o  noon  o  o  o  o  noon  on o non 


SUBROUTINE  KUWSRT (A,B,k,N) 


*  SUBROU I  INF  kUrtSRT  SORTS  ROWS  OF  A  MATRIX  A  ACCORDING  TO  DESCENDING 
ORDER  OF  MATRIX  B.  STORE  S0R1ED  A  INTO  R. 


DIMENSION  A1810),  B(blO),  R(lb20) 

COMMON  / SRT  COM/  1 ,  1 2 , I  SORT , KS A VE , SA VER , I R , 1  A 
LEVEL  2,  A,6,R, I,N, 12, I  SORT , RSAVE , SAVER , IK, IA 

*  MOVt  SORTING  KEY  VECTOR  TO  FIRST  COLUMN  OF  OUTPUT  MATRIX  AND 

BUILD  ORIGINAL  SEQUENCE  LIST  IN  SECOND  COLUMN 

DU  1000  1  =  1  ,  N 
R  ( 1 1  =  B  ( I  ) 

12  =  1  +  N 
R  ( 1 2)  =  I 
1000  CONTINUE 

*  SORT  ELEMENTS  IN  SORTING  KEY  VECTOR  (ORIGINAL  SEQUENCE  LIST  IS 

RESEQUENCED  ACCORDINGLY) 

2000  CONTINUE 
ISORT  =  0 
DO  4  000  1  =  2 f  N 

IF  ( R ( I )  -  R(I-l))  4000,  4000,  3000 
3000  CONTINUE 

ISORT  =  ISORT  +  1 
RSAVE  =  R  ( 1 ) 

R ( I )  =  R(I-l) 

R ( I - 1 )  =  RSAVE 
12  =  1  +  N 
SAVER  =  R ( 12) 

R ( 1 2 )  =  R ( 1 2- 1  ) 

R ( 1 2- 1 )  =  SAVER 
4000  CONTINUE 

IF  (ISORT)  2000,  5000,  2000 

*  MOVE  ROWS  FROM  MATRIX  A  TO  MATRIX  R  (NUMBER  IN  SECOND  COLUMN  OF  R 

REPRESENTS  NOrt  NUMBER  OF  MATRIX  A  fO  BE  MOVED) 

5000  CONTINUE 

DO  6000  I  =  1,N 

*  MOVE  ELEMENT  TO  OUTPUT  MATRIX 

12  =  I  +  N 
R ( I )  =  A (  R ( 12)  ) 

8000  .CUNT  I NUE 
RETURN 
END 
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SUBKUUT 1 N  E  SCRIB2 


- - 

c 

C  *  SCRIBE  PLOTS  GRAPHS  FOR  MUZZLEDE VICE  DA  1  A  IN  PANS 
C  USING  CALCOMP  PLOTUNb  PACK AGt 

C 

- - 

COMMON  /  SCR1B2  /  XARRAY  (100)  ,  YARRAY  (100)  ,  DX1  ,  NDPTS  ,  IMPP1 

REAL  DX  ,  XSCALE  ,  YSCALE  ,XLNGTH 

REAL  YLNGTH  ,  XOR  ,  YOR  ,  XBASE  ,  YBASE  ,  FACTOR 

INTEGER  I  SENT (10)  ,  L AbEL  (4 )  ,  1UNIT 
COMMON  /  LABEL  /  KTB(2«)  ,  JSTART  ,  I 
LEVEL  2,  KFh  ,  TSTaRT  ,  I 
C 

DATA  LABEL/ ' BOOT S ', '6390 X39B3 ', 'M392  DATA'/ 

C 

C  TEST  ARGUMENTS  FUR  ERRORS 

C 

IF  (NDPTS  .GT.  100)  CALL  HERROk  (93) 

C 

C  DEFINE  VARIABLES  USED  BY  PLOTTING  ROUTINES 

C 

XLNGTH=11  . 

YLNGTH=8.5 
XS12E=B.O 
YSIZE=6.0 
XBASE=2 . 0 
YBASE=1 .0 

CALL  FIXSCA  (XARRAY 11), NDPTS, XSIZE, X  SC ALE, XOR, XMAX,DX) 

CALL  FIXSCA ( YARRAY (1 ) , NDPTS, YSIZE, YSCALE, YOR, Y  MA  X , D Y ) 

XLNGTH=1 1  . 

YLNGTH=8.b 
XS1ZE=8.0 
YSIZE=6.0 
XBASE=2.0 
YBASE= 1 .0 

CALL  FIXSCA (XARR AY (1 ) , NDPTS, XSIZE, X SC ALE, XOR, XMAX,DX) 

CALL  FIXSCA ( YARRAY (1), NDPTS, YSIZE, YSCALE, Y UR, YM AX, DY) 

FACTOR  =1.  • 

IUNIT  =  10 
C 

C  BRL  PLOTTING  ROUTINES  (INITIAL) 

C 

CALL  PLTBEG ( XLNGTH, YLNGTH , F AC T OR , I  UN  I T , L ABEL ) 

CALL  PLT SC A(XBASE,YBASE, XOR, YOR, XSCALE, YSCALE) 

C 

C  PLOT  THE  CURVES 

C 

CALL  PL TDTS(1,0, XARRAY, YARR AY, NDPTS,0) 

C 

C  PLOT  AXES  AND  LABELS 

C 

MODE  =  4 
ONE  =  1. 

CALL  PLTAXS(OX,DY, XOR, XMAX, YOR,  YMAX,  MODE) 

CALL  LABELA (DX , DY , XOR, XMAX, YOR,  YMAX  ,  ONE,  ONE) 

C 
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YMAXlsYMAX+1 .  0*YSCALE 
YMAX 1=YMAX+1 .0* YSCALE 
ENCODE  11UU,  130,  1  SEN  1 ) 
130  EUWrtAT  (IX,  P4I2,  •>') 
CALL  PL1SYM(0.15,ISENT,0 
CALL  PLTSYM(0. lb, ISENT  ,0 

CALL  PLTPGE 

RETURN 

END 
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. , XUR, YMAX1  ) 
. , XOK, YMAX1  ) 
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ooo  ooo  ooo  noon  o  o  o  o  o  o 


OVERLAY  (  J  ,  b  » 3  ) 
PRUGRAM  SLR1BL 


*  SCR  1 BL  PLOTS  UP  TO  THREE  OVERLAPPING  GRAPHS  FOR 

MIJZ ZLEDE V ICE  DATA  USING  CALCOMP  PLOTTJNG  PACKAGE 


COMMON  /  SCRIB  /  XARRAY  15000)  ,  YARRY1  15000)  ,  YARRY2(bO00) 
A  YARRY3(5000)  ,  NARRAY  ,  0X1  ,  NOPTS  ,  NPPI 

LOGICAL  TWICE 

REAL  DX  ,  XSCALE  ,  YSCALE  , XLNGT  H 

REAL  YLNGTH  ,  XOR  ,  YOk  ,  XuASE  ,  YBASE  ,  FACTOR 

INTEGER  ISENT(IO)  ,  LAbEL(<4)  ,  IUNIT 
COMMON  /  LABEL  /  KTB  (2*4 )  ,  T  ST  ART  ,  I 

LEVEL  2,  KTB,  TSTART,  I 

DATA  LABEL/ *  BOOTS'  ,  'B390  X3963 ', 'M392  DATA1/ 

TEST  ARGUMENTS  FOR  ERRORS 
IF  (NOPTS  .GT.  5000)  CALL  PERROR  (84) 

IF  (  (NARRAY  .LT.  1)  .OR.  (NARRAY  .GT.  3)  )  CALL  PERROR  (83) 
OEF1NE  VARIABLES  USED  BY  PLOTTING  ROUTINES 


XLNGTH=1  1  . 

YLNGTH=8  .  5 
XSIZE=8.0 
Y  S I ZE=6 . 0 
XBASE=2.0 
YBASE  =  1 .0 

CALL  FIXSCA(XARRAY(1),  NOPT S , X S I L E , X SC ALE , XOR , XMAX,DX) 
CALL  FIXSCA (YARKY1 (1) , NOPTS, Y S I Z E , Y SC  ALE , Y OR , Y M A  X , D Y ) 
XLNGTH= 1 1 . 

YLNGTH=6.5 
XS1ZE=8.0 
YSIZE=6.0 
XbASE=2.0 
YBASE  =  1 .0 

CALL  FIXSCA (XARRAY (1)  ,  NOPTS, XSIZE,XSCALE,XOR,XMAX,DX) 
CALL  FIXSCA ( YARRY1 (1 ) , NOPTS, YS1ZE, YSCALE, YOR, YM AX, OY) 
FACTOR  =  1. 

IUNIT  =  10 
TWICE  =  .FALSE. 

BRL  PLOTTING  ROUTINES  (INITIAL) 

CALL  PLT BEG ( XLNGT H,  YLNGTH,  FACT  OR, IUNIT, LABEL) 

CALL  PLTSCA (XBASE, YBASE, XUR, YOR, XSCALE, YSCALE) 

PLOT  THE  CURVES 

GO  TO  (3000,2000,1000)  ,  NARRAY 

1000  CONTINUE 

CALL  PL TDTS( 1,0,  XARRAY, YARRY3, NOPTS, 0) 

2000  CONTINUE 
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LALL  PL  TUTS  (4,0,  X ARRAY  ,  Y  ARR  Y  2 ,  NOP  T  S ,  0  ) 
3000  CUNT INUt 

CALL  PL  I  UTS (1,0, X ARRAY, YARRY 1 , NDPTS, 0) 
PLOT  AXES  AND  LABELS 


MODE  =  4 
ONE  =  1. 

CALL  PLTAXStUX, DY, XOk, XMAX, YON, YMAX , MODE) 
CALL  L  ABEL  A  (f)X,OY,  XOk,  XMAX,  YON,  YMAX,  ONE,  ONE) 
C 

YMAX1=YMAX+1.0*YSCALE 
YMAX1=YMAX+1 .0*YSCALE 
ENCODt  (100,  130,  1SENT)  KTB 

130  FORMAT  (IX,  2412,  '>') 

CALL  PL TSYM(0. 15,1  SENT, 0., XOk, YMAX1) 

CALL  PLTSYM(0.15,lSENT,0.,X0k,YMAXl) 

C 

C  IF  (TwICE)  GO  TO  500 

CALL  PLTPGE 
C  500  CONTINUE 

TWICE  =  .TRUE. 

RETURN  * 

END 
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UVfckLAY (G,b, 0) 

PROGRAM  SMOOTH 

C - - - - 

C 

C  *  SUBROUTINE  SMOOTH  ALLOTS  THE  USER  TO  SELECT  A 
C  SMOOTHING  OPTION. 

C 

C - FALCON  R&D  MEK  7-  9-79 

COMMON  /  CONTRL  /  LUIN  ,  LUOUT  ,  LUDATA  ,  LURUN  ,  LUREPT 
COMMON/C ONTRL/IFSVRP, ONLINE, I ANOPT , OPTION, IFBACH 
CUMMON  /  CONTRL  /  USRNAM  ,  PURPOS 

COMMON  /  SSOATA  /  NSSPTS,  S0T1ME  ,  SUBSTR(5000)  ,  SSTIMEC5000) 
INTEGER  ONLINE  ,  OPTION  ,  USkNAM(2)  ,  PURP0S(6) 

INTEGEk  SMOPT 

LEVEL  2,  NSSPTS,  SOT  I  ME ,  SUBSTk,  SSTIME 

LEVEL  2,  L 0 IN, LUOUT, LUDATA, LUk UN, LUREPT , I FS VRP, ONL I NE , I ANOPT 
LEVEL  2,  OPTION, USRNAM, PUkPOS, SMOPT, NSMO 
LEVEL  2, IFBACH 

COMMON  /  COMSMO  /  SMOPT  ,  NSMO 
C 

C  *  USER  SELECTS  A  SMOOTHING  OPTION  (BATCH) 

C 

READ  (LUIN,  9002)  SMOPT 
*  SELECT  NUMBER  OF  POINTS  FOR  SMOOTHING 


2000  CONTINUE 

WRITE  (b,  9003)  SMOPT 
IF  (SMOPT  .EU.  0)  GO  TO  9000 

IF  (SMOPT  .LI.  1  -OR.  SMOPT  .GT.  2)  CALL  PERR0R(99) 
GO  TO  (3000,  9000),  SMOPT 
3000  CONTINUE 

*  11-POINT  MEDIAN  -  SMOOTHING  OF  SUBSTk 

CALL  MED  1 1 
GO  TO  6000 
9000  CONTINUE 
C 

C  *  NSMO-POINT  AVERAGE  -  SMOOTHING  UF  SUBSTR 
C 

NSMO  =  11 

CALL  MOVAVE  (SUBSTR,  NSSPTS,  NSMO) 

GO  TO  8000 
8000  CONTINUE 
9000  RETURN 

9002  FORMAT  (12) 

9003  FORMAT  (IX,  20HSM00THING  OPTION  IS  ,  12) 

END 
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SUrtKOUT INt  MED  11 


*  SUBROUTINE  MtUll 

11-POINT  MEDIAN  -  SMOOTHING  OF  SUBSTR 

_ _ _ _ _ -  —  -FALCON  R&D  NLA  7-  9-79 

COMMON  /  SSDATA  /  NSSPIS  ,  SOT  1 ME  ,  SUBSTRC5000)  *  SSTIMEC5000) 
COMMON  /  CUMMU  /  L ,  Z,  ZM,  1U 
DIMENSION  Zlll) 

COMMON  /  GENSTR  /  XC5000),  LL,  1,  LS,  LN,  K,  KS 
LEVEL  2,  NSSPTS ,  SOTIME,  SUBSTR,  SS1IME,  X 
LEVEL  2,  LL,  L,  I,  L,  ZM,  IU,  LS,  LN,  K,  KS 
L  =  11 
LL  =  ,5  *  L 
DO  1000  I  =  1 , L 
ZCI)  =  SUbSTRCI) 

1000  CONTINUE 

CALL  MED  I  AN ( Z , L , ZM) 

00  1500  I  =  1 , LL 
X(I)  =  ZM 
1500  CONTINUE 

00  2000  I  =  1 , L 

IU  =  NSSPTS  -  I  +  1 
ZCI)  -  SUBSTRCI) 

2000  CONTINUE 

CALL  MEDIANCZ,L,  7.M) 

00  2500  1  =  1 , LL 

IU  =  NSSPTS  -  I  +  1 
XCIU)  =  ZM 
2500  CONTINUE 

LS  =  LL  +  1 
LN  =  NSSPTS  -  LL 
00  3000  K  -  L S , L N 
KS  =  K  -  LL  -  1 
00  4000  I  =  1,L 

Z(D  =  SUBSTRCKS  +1) 

4000  CONTINUE 

CALL  ME0IANCZ,L,ZM) 

X  CK)  =  ZM 
3000  CONTINUE 

DO  5000  1  =  1, NSSPTS 
SUBSTRCI)  =  XU) 

5000  CONTINUE 
RETURN 
END 
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SUBROUTINE  MEDIAn(Z,N, ZM) 


c 


*  SUBROUTINE  MEDIAN  RETURNS  THE  MEDIAN  ZM  OF  N  NUMBERS 
IN  ARRAY  Z 


- - - - FALCUN  RND  NEA  7-  9-79 

DIMENSION  Z ( N ) 

COMMON  /  COMMED  /  L,  J,  AMIN,  IMIN,  I 
LEVEL  2,  Z,  N,  ZM,  L,  J,  AMIN, IMIN, 1 

l  =  n  /  a  +  1 

DU  1  000  J  =  1  ,  L 
AMIN  =  10000000. 

DO  20  00  I  =  1  ,  N 

IF  (Z(I)  .GE.  AMIN)  GO  TO  2000 
IMIN  =  1 
AMIN  =  Z  C 1  ) 

2000  CONTINUE 

IF  (J  .fcy.  L)  GO  TO  3000 
ZC1MIN)  =  1000000. 

1000  CONTINUE 
3000  CONTINUE 

ZM  =  Z(IMIN) 

RETURN 

END 
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SUtiROUT  i  NE  MUVAVE  (S,  Nr  L) 

C - * - 

C 

L  *  SUrtKUUTINE  MOVAVE 

C  L-POINT  AVERAGE  -  SMOOTHING  OF  ARRAY  S  WITH  N  ELEMENTS 

C 

C - FALLON  k&D  NEA  7-  9-79 

DIMENSION  S ( N  7 

COMMON  /  COMMOV  /  SS,  SL,  IL,  SA,  SC 

COMMON  /  GENSTR  /  X (5000) ,  LLr  lr  LS,  LN,  K,  KS 

LEVEL  £  ,  S  r  Nr  Lr  X 

LEVEL  £ ,  LL,SS,SL, IL, 1 , SA , SC , LS , LN , K , KS 
LL=  .b  *  L 
SS  =  0. 

SL  =  0. 

DO  1O00  I  =  1 , LL 
IL  =  N  -  I  +  1 
SS  =  SS  +  S(I) 

SL  =  SL  +  SUL) 

1000  CONTINUE 

SA  =  SS  /  FLOAT (LL) 

SC  =  SL  /  FLOAT (LL) 

DO  £000  I  =  1,LL 
IL  =  N  -  ]  t  1 
X ( I )  =  SA 
X  (  I  L  )  =  SC 
£ 0 0 0  CONTINUE 

LS  =  LL  +  1 

LN  =  N  -  LL 

DO  9  0  0  0  K  =  L S , L N 

X ( K )  =  0. 

KS  =  K  -  LL  -  1 
DO  3000  I  =  1  ,  L 

X (K)  =  X (K)  +  S(KS  +  I) 

3000  CONTINUE 

X ( K )  =  X ( K )  /  FLOAT (L) 

9000  CONTINUE 

DO  5000  K  =  1 , N 
S  (  K  )  =  X  ( K  ) 

5000  CONTINUE 
RETURN 
END 
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OVERLAY  (L,  7,  0) 

PROGRAM  L/LKIV 

C - - - - - - 

C 

C  *  SUBROUTINE  UEKIV  CALCULATES  THE  FIRST  ANl)  SECOND 
DERIVATIVES  OF  ENTIRE  DATA  STkING.  PLOTS  DAIA 
STRING  AND  DERIVATIVES. 


COMMON  /  SDATA  /  NS1PIS  ,  DT 1  ME  ,  STR1NG(5000)  ,  STIME(5000) 

COMMON  /  UERIVS  /  STRDOT(5000)  ,  STRDD(SOOO)  ,  PLOTDP 

COMMON  /  SCR  1 64  /  XARKAY(SOOO)  ,  YARRY1 (5000)  ,  YAKRY2(5000)  , 

*  YARRY3(bOOO)  ,  I)X1  ,  NDPTS  ,  NPPI 

COMMON  !  LAbEL  /  KT6(24)  ,  TSTAR1  ,  I 

COMMON  /  CLASCM  /  ICLASS  ,  NBEGBL 
COMMON  /  COMDER  /  NS  T  M2  ,  CD  ,  DT  ,  CDO 

LEVEL  2,  NSTPTS,  DTIME,  SIRING,  ST1ME,  STRDDT,  STRDD,  KTb,TSTART 

*  ,  I,  PLOT OP,  ICLASS, NuEGBL,NSTM2, CD, DT, COD 
INTEGER  PLOTDP 

*  WRITE  PAGE  HEADING  FOR  DERIVATIVES 

wR I TE ( b , 90  04 ) 

*  PLOTOP  =  PLOT  OPTION  FOR  DERIVATIVES 

=  1  ON 
=  0  OFF 

RE  A  D  ( 5 , 90  0  2 )  PLOTOP 
WRI TE  (6, 90031  PLOTOP 
DO  1000  I  =  1, NSTPTS 
STRDtJ T  (I)  =  0. 

STRDD(I)  =  0. 

1000  CONTINUE 

DT  =  DTIME  *  .001 
CD  =  1.  /  (12.  *  DT) 

COD  =  CD  /  DT 

NS T M2  =  NSTPTS  -  2 

*  ICLASS  INDICATES  SATURATION  LEVEL: 

ICLASS  =  -1,  HEAVY  SATURATION,  EPOCH  1  DETECTION 
AND  REGIME  1  ANALYSIS  ONLY. 

ICLASS  =  0,  NO  SATURATION,  FULL  ANALYSIS. 

ICLASS  =  1,  LIMITED  SATURATION  (CLIPPING), 

FULL  ANALYSIS. 

IF  (ICLASS  .EQ.  -1)  NSTM2  =  NBEGBL  -  2 
DO  2000  I  =  3 , NSTM2 

STRDOT (1)  =  CD  *  (6.  *  (STR1NG(1+1)  -  STRING! 1-1 ) ) 

*  -  STRING! 1+2)  +  STRING  (1-2)) 

STRDD(I)  =  CDD  *  (16.  *  (STRINGU  +  l)  +  STRlNG(I-l)) 

*  -  30.  *  STRING (1)  -  STRING( 1+2)  -  STRING(I-2)) 

2000  CONTINUE 

*****  ****************************** 
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C  *  wkITE  (JUT  FIRST  AND  SECOND  DERIVATIVES 
C 

rtRlTE(b,9005) 

*.RITE(b,900b) 

DU  250  0  1  =  1,NSTP1S 

nKITE(b,9001)  I,STIME(I) , ST  RING ( 1 ) , STKUOT ( I ) , STkDD ( 1 ) 

S2SOO  CONTINUE 

C*********************************** 

c 

IF  (PLOTOP  .  EQ  .  0)  60  TO  <(000 
C 

DO  3000  1  =  1 , NST PT S 
XAkRAY(l)  =  STIME(I) 

YARKYl(I)  =  STRING(I) 

Y  ARk Y2  ( I )  =  STRDOT (I) 

$  Y  Akk  Y3 ( I )  =  STRDO(I) 

3000  CONTINUE 

DX1  =  DTIME 
NDPTS  =  NSTPTS 
NPP I  =  80 

CALL  OVERLAY (1HM, 7,1) 

<4000  CONTINUE 

C*********************************** 
90  0  1  FORMAT  (19X,1<4,2X,F9.5,2X,F9.5,2X,E13.7,2X,E13.7) 

9002  FORMAT (12) 

9003  FORMAT  (IX,///,  10X,2<4HPL0T  OPTION  FOR  UERIV  IS,5X,12) 

90  0 <4  FORMAT  (lHl,///,20X,54(lh*),/,20X,lH*,S2X,lH*,/,20X,lH*,15X, 

*  22H DERIVATIVE  C ALCUL AT  I  ON, 1 5 X , 1 H* , / , 2 0 X , 1 H* , 52X , 

*  1H*,/,20X,54(1H*)) 

9005  F0RMA1 (1X,/,10X,29HFIRST  AND  SECOND  DERIVATIVES  , 

*  3 1  HARE  CALCULATED  BUT  NOT  PRINTED.) 

90  0b  FORMAT  Cl  X,b(/)  ,  1  8  X  ,  ShPO  I  NT  ,  <4  X  ,  4  HT  IME  ,  7  X  ,  4HD  A  T  A  ,  8  X  , 

*  5HFIRST ,9X,bHSECUN0,//) 

RETURN 

END 
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0VERLAy(M,7,1) 
PROGRAM  SCK1H4 


*  SCR1B4  PLOTS  THREE  GRAPHS  OF  THt  ORIGINAL  DATA  STRING  AND  ITS 
FIRST  AND  StLUNU  DfcRl  VAT  I VF.  S  USING  CALCOMP  PLOTTING  PACKAGE 


COMMON  /  SCRIBE  /  XARRAY (5000)  ,  YARRY 1  (5000 1  ,  YARRY2(5000)  , 

A  YARRY3(5000)  ,  DX1  ,  NDPTS  ,  NPPI 

REAL  DX1  ,  XSCALE  ,  YSCALE  ,  XLNGTH 
REAL  YLNGTH  ,  XOR  ,  YOR  ,  XBASE  ,  YBASE  ,  FACTOR 
INTEGER  lSENT(lO)  ,  LABEL  (4  )  ,  IUNIT 

COMMON  /  LABEL  /  KTb ( 24)  ,  TSTART  ,  I 

LEVEL  2,  KTB  ,  TSTART  ,  I 

DATA  LABEL/ 'BOOT S B390 ' , ' X3983 ' , ' M392  DATA'/ 


IF  (NDPTS  .GT.  5000)  CALL  PERR0R(95) 

DEFINE  VARIABLES  USED  TO  PLOT  YARRY1 


MODE  =  0 
ONE  =  1. 


XLNGThr 1  l  . 

YLNGTH =8. 5 

XS1ZE=S.0 

YSIZE1=4. 

YSI ZE2=  1  .5 
XBASE=2. 

YBASE=3.0 

CALL  FIXSCA (XARRAY (1) , NDPTS, XSIZE, XSCALE, XOR, XMAX,DX) 
CALL  FIXSCA (YARRY 1 (1 ) , NDPTS, YS I  Zb 1 , Y SC ALE , Y OR , YM AX , D Y ) 
XLNGTH=1 1 . 

YLNGTH=B.5 

XSIZE=8.0 

YSIZE1=4. 

YSIZE2-1 .5 
XBASE=2 . 

YBASE=3.0 

CALL  FIXSCA (XARRAY (1) , NDPTS, XSIZE, XSCALE, XOR, XM A X,DX) 
CALL  FIXSCA (YARRYl(l) , NDPTS, YSIZE1, YSCALE, YOR, YMAX,DY) 
FACTOR  =  1. 

IUNIT  =  10 


PLOT  DATA  STRING 


C 


CALL  PLTBEG(XLNGTH, YLNGTH, FACTOR, IUNIT, LABEL) 
CALL  PLT SC A (XBASE, YBASE, XOR, YOR, XSCALE, YSCALE) 
CALL  PLTOTS( 1,0, XARRAY, YARRY 1 , NDPTS, 0) 

CALL  PLTAXS(DX,DY,XOR,XMAX,YOR, Y MAX, MODE) 

CALL  LABELA(0. ,DY, XOR, XM A X , Y OR , YM A  X , ONE , ONE ) 
CALL  L ABEL A (0. , DY, XOR, XMA X , Y OK , Y M A X , ONE , ONE ) 
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C  UEHNE  VARIAbLES  USED  TO  PLOT  Y  ARR  Y2 

C 

CALL  F  IXSCA  (  YARRY2U  )  ,NUPTS,  YS1ZE2,  YSC  ALE  ,  YOk ,  Y  MA  X  ,  D  Y  ) 
YbASE=l  . 

CALL  F IXSCA(YARRY2(1 ) ,NDP1S, YSIZt2, Y SCALE#  YOR, YMAX , UY ) 

Y  ti  A  S  E  =  1  . 

PLOT  FIRST  DERIVATIVE 

CALL  PL T SC A (XBASE, Y BASE, X OR, Y OR, X SC ALE, YSC ALE) 

CALL  PLTDTSU,  0,X  ARRAY,  Y AKRY2, NDPT S,  0) 

CALL  PLT  AXS (DX, DY, XOk, XMAX, YOR, YMAX, MODE) 

CALL  PLT AXS (OX , DY , XOR , XMA X , YOR , YMAX , MODE ) 

CALL  LABEL A ( DX , D Y , XOK , XM A X , Y OR , Y M AX , ONE , ONE ) 

DEFINE  VARIABLES  USED  TO  PLOT  YARRY3 

CALL  F1XSCA(YARRY3(1) , NDPTS,b., Y $3, SMI  NOD, SMAXDD, YELD3) 
S>  Y  SC  ALE  =  YS3 

5>  DY  =  YELD3 

$  YOR'  =  SMINDD 

$  YMAX  =  SMAXDD 

3>  YBASE  =  22. 

$  XMAX  =  XOR 

C 

C  PLOT  SECOND  DERIVATIVE 

C 

$  CALL  PLT SCA (XBASE, YBASE, XOR, YOR, XSCALE, YSCALE ) 

S  CALL  PLTDTSU, 0,XARRAY,  Y  ARR  Y  3 ,  NDP1  S ,  0  ) 

5  CALL  PLTAXS(0.,DY, XOK, XMAX, YOR, YMAX, MODE) 

3>  CALL  LABEL  A  (DX,  DY  ,  XOR,  XMAX,  YOR,  YMAX,  ONE,  ONE) 

C 

YMlNl=YOR+(YSCALE*b.3) 

YMlNl=Y0R+(YSCALE*b.3) 

ENCODE (100,  130,  1  SENT )  KTB 

130  FORMAT (IX, 2412, ' >  '  ) 

CALL  PLTSYM(0. lb, I  SEN T , 0 . , XOR , YM I N 1 ) 

CALL  PLT5YM(0.1b,ISENT,0.,XOR,YMlNl) 

CALL  PLTPGE 

RETURN 

END 
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OVERLAY (N, 1 0, 0) 

program  detect 


c - -- - - - --- - 

c 

C  *  SUBROUTINE  DETECT  TRIES  TO  LOCATE  THE  FIRST  EPOCH 
C  AND  THE  SECOND  EPOCH. 

C 

- - 

COMMON  /  SO AT  A  /  NSTPTS  ,  DTIME  ,  STRING(5000J  ,  STIMEC5000) 

COMMON  /  DERIVS  /  STRDOT ( 5000 1  ,  STRDD(5000)  ,  PLOTOP 

COMMON  /  COMDET  /  A  ( b  1  0 )  ,  B  ( 8 1  0  )  ,  R  C 1 620)  ,  NROW ,  1 ,  IMMAXOV ,  STRMA  X ,  IP 

*  ,  N(200),  STR ( 200 ) ,  STRM1N  f  b  ACK  ,  CURR 
COMMON  /  DTECT  /  EPUCH1  ,  EP0CH2 

COMMON  /  CLASCM  /  1CLASS  ,  NBEGBL 
INTEGER  EPOCH  1 ,  EP0CH2  ,  PLOTOP 

LEVEL  2  ,  A,  8,  R,  NROw,  NSTPTS,  DTIME,  STRING,  STIME,  STRDOT, 

*  STROD,  PLOTOP,  I,  NMAXDV,  STRMAX,  IP, 

*  n, STR, STRMIN, EPOCH1 , EP0CH2, I  CL  ASS, NBEGBL 
LEVEL  2, BACK, CURR 

C 

C  *  WRITE  PAGE  HEADING  FOR  DETECT 
C 

WRITE (6, 9003) 

C 

C  *  FIND  THE  MAXIMUM  OF  STRDOT 
C 

STRMAX  =  -99999. 

DO  2000  I  =  1, NSTPTS 

IF  (ST  ROOT ( I )  -  STRMAX  .LE.  0)  GO  TO  1000 
STRMAX  r  STRDOT (I) 

NMAXDV  =  I 
1000  CONTINUE 
2000  CONTINUE 

C 

C  *  SORT  STRING 
C 

NROrt  =  200 

DU  3000  I  =  1 ,  N  R  0  W 

Adi  =  float  (  nmaxdv  -  (nrow  -  i)i 

B ( I )  =  STRING (NMAXDV  -  (NROrt  -  I)) 

3000  CONTINUE 

CALL  ROWSRT  (A,8,R,NR0w) 

C 

C  *  FIND  THE  EPOCH 
C 

C  *  IN  AN  INTERVAL  OF  200  POINTS  PRECEDING  MAX  DERIV, 

WORN  BACKWARDS  TO  FIND  LAST  POINT  GREATER  THAN 
POINTS  PRECEDING  IT. 

DO  9000  I  =  1 , NROW 
IP  =  NMAXDV  -  I  +  1 

IF  ( ABS (R ( I )  -  FLOAT ( IP) )  .GT.  9.)  GO  TO  5000 
9000  CONTINUE 
5000  CONTINUE 

EPOCH  1  =  IP  +  1 
C 
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C  *  WORK  BACKWARDS  TO  FIND  LASl  muinu  i  uw  x  i«ll  i 
C  OECkEASIwb  POINT.  THAI  IS  EP0CH1 . 

C 

CURR  =  EPOCH  1 
5400  CONTINUE 

BACK  =  CuRW  -  1 

IF  (BACK  .LE.  1.)  CALL  PERR0k(65) 

IF  (STRING(BACK)  .GE.  STRING ( CURR ) I  GO  TO  5800 
CURR  =  BACK 
GO  TO  5400 
5800  CONTINUE 
C 

EP0CH1  =  CURR 

C  **********************************************************************' 


c 

c 

c 

c 

c 

c 

c 

c 


*  ICLASS  INDICATES  SATURATION  LEVEL: 

1CLASS  =  -1,  HEAVY  SATURATION,  EPUCH  1  DETECTION 
AND  REGIME  1  ANALYSIS  ONLY. 

ICLASS  =  0,  NO  SATURATION,  FULL  ANALYSIS. 

ICLASS  =  1,  LIMITED  SATURATION  (CLIPPING), 

FULL  ANALYSIS. 

IF  (ICLASS  .EQ.  -1)  GO  TO  9000 


*  IDENTIFY  FIRST  200  VALUES  AFTER  NMAXDV  OF  STRING 

NROW  =  200 

DO  6000  I  =  1 ,  NROW 

N ( I )  =  NMAXDV  +  (I  -  1) 

STR(l)  =  STRING (N(I) ) 

6000  CONTINUE 

*  FIND  FIRST  MINIMUM  OF  STR(I) 

6500  CONTINUE 

STRMIN  =  9999999. 

IP  =  0 

DO  8000  I  =  1 , NROW 

IF  (STR(I)  .GE.  STRMIN)  GO  TO  7000 
STRMIN  =  STR(I) 

IP  =  N  (  1 ) 

7000  CONTINUE 
8000  CONTINUE 

EP0CH2  =  IP 


*  CHECK  IF  3  PRECEDING  POINTS  TO  EPOCH2  ARE 
MONOTONICALLY  DECREASING. 

CURR  =  EP0CH2 
DU  8300  1  =  1,3 
BACK  =  CURR  -  1 

IF  (STRING(BACK)  .LE.  STRING (CURR) )  GO  TO  8500 
CURR  =  BACK 
8300  CONTINUE 
C 

C  *  EP0CH2  IS  VALID 


PROGRAM  LISTING  (Continued) 
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onoo  o  ooo 


c 

GU  TD  d  7  0  0 
6500  CONTINUE 


*  EP0CH2  15  INVALID,  REPEAT  PROCESS  10  THAT  POINT 


NRUW  =  CURR  -  NMAXDV 
GO  TO  6500 
6700  CONTINUE 


WKITE(6,9001)  EPOCH  1, ST  1  ME ( EPOCH  1) , ST  RING (EPOCH  1 ) 
rtRITE (6,9002)  EPOCHS, ST  1  ME (EPOCHS) , STRING (EPOCHS) 


*  FIND  LOCAL  EXTREMA  IN  DERIVATIVE 

NST=NMAXDV-40 
IF (NST.L1 .2)NST=2 
NSP=NST+120 
CR=STkMAX* . 06 
EXTR=STRDOT (NST ) 
rtRITE (6, 9006) 

DO  6900  1=NST,NSP 
Al  =  STRDOT 1 1 - 1 ) 

A  2  =  S  T  R  D  0  T  (I) 

A3  =  STR0UT ( I  +  l ) 

lF(Al.LE.A2.AND.AS.GE.A3.AND.AbSlA2-tXTR).GE.CR)GU  TO  6600 
IF ( A  1 .GE. AS. AND. A2.LE. A3. AND. AdS ( A2-EXTR) .GE.CR)GO  TO  6650 
GO  TO  8900 

6600  WKITE(6,9004)ST1ME(1), AS 
EXTR=A2* 

GO  TO  8900 

8650  rtRlTE(6,9005)STlME(I),A2 
EXTR=A2 
8900  CONTINUE 

9001  F'URMAT  ( 1  X, ////,  10X,  18HFIRST  EPOCH  IS  AT  ,15,3X, 

*  10HWITH  TIME  ,F10.4,5H  AND  , 

*  llhwITH  VALUE  ,F10.4) 

90 OS  FORMAT ( 1 X, /, 10X, 18HSECOND  EPOCH  IS  AT  ,I5,3X, 

*  10HWITH  TIME  ,F10.4,5H  AND  , 

*  11HWITH  VALUE  ,F10.4) 

9003  FORMAT (1H1,///,20X,54(1H*),/,20X,1H*,52X,1H*,/,20X,1H*,17X, 

*  1 9HDETECTIUN  OF  EPOCHS , 1 6X , 1 H* ,/, SO X , 1 H* , 52X , 1 H* , 

*  / , SO  X , 54 ( 1 H* ) ) 

9004  FORMAT (20X,F 10.4, 1 0 X , F 1 1 .0) 

9005  FORMAT (20X,F10.4,30X,F11.0) 

9006  FORMAT (IX,////, SOX, 54(1H*),/, SOX, 1H*,52X,1H*,/,20X,1H*,13X, 

*  27HLOCAL  EXTREMA  IN  DER I V AT  I VE , 1 SX , 1 H* , / , SOX , 1 H* , 5SX , 1 H* , / , 

*  20X,54(1H*),////,25X,4HT1ME, 14X, 7HMAXIMUM, 13X, 7HM1NIMUM, /) 
9000  RETURN 

ENO 
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SUBROUTINE  KUASk 1  l  A, B, R, N) 

- - 

C  ■ 

C  *  SUBROUTINE  RU'ASRf  SORTS  ROWS  OF  A  MATRIX  A  ACCORDING  TO  DESCENOINt 
C  UNDER  OF  MATRIX  B.  STOKE  SORTED  A  INTO  R. 

C 

- - 

DIMENSION  A(810),  B(810),  R(lbEO) 

COMMON  /SRTCOM/  I ,  I'd,  ISORT,  RSAVE,  SAVER,  IR,  IA 
LEVEL  2,  A, B,R, I , N, 12, ISORT ,RSAVE, SAVER. IR, I A 
C 

C  *  MOVE  SORTING  KEY  VECTOR  TO  FIRST  COLUMN  OF  OUTPUT  MATRIX  AND 
C  BUILD  ORIGINAL  SEQUENCE  LIST  IN  SECOND  COLUMN 

C 

DU  1000  I  =  1,N 
R  (I )  =  Bill 
IE  =  1  +  N 

RUE)  =  1 
1000  CONTINUE 
C 

C  *  SORT  ELEMENTS  IN  SORTING  KEY  VECTOR  (.ORIGINAL  SEUUENCE  LIST  IS 
C  RESEUUENCED  ACCORDINGLY) 

C 

E000  CONTINUE 
ISORT  =  0 
DO  4  0  0  0  1  =  2 , N 

IF  (RID  -  R  ( 1  - 1 ) )  4000,  4000,  3000 
3000  CONTINUE 

ISORT  =  ISORT  +  1 
RSAVE  =  k(I) 

R ( I )  =  R(l-l) 

R(I-l)  =  RSAVE 
IE  =  1  +  N 
SAVER  =  RUE) 

RUE)  =  RUE-1) 

RUE-1)  =  SAVER 
4000  CUNTINUE 

IF  (ISORT)  E000 ,  5000,  EOOO 
C 

C  *  MOVE  kOAS  FRUM  MATRIX  A  TO  MATRIX  R  (NUMBER  IN  SECOND  COLUMN  OF  k 
C  REPRESENTS  ROW  NUMBER  OF  MATRIX  A  TO  BE  MOVED) 

5000  CONTINUE 

DO  80  00  I  =  1  ,  N 
C 

C  *  MOVE  ELEMENT  TO  OUTPUT  MATRIX 
C 

IE.  =  1  +  N 

R  ( I  )  =  A  (  RUE)  ) 

RUE)  =  B  (  RUE)  ) 

8000  CONTINUE 
RETURN 
END 
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o  o  o 


OVERLAY  (P, 1 1 , 0) 

PROGRAM  PEG2AN 

c  — - - - 

C 

C  *  SUBROUTINE  REG2 AN  --SPECIAL  ANALYSIS  UF  REGIME  II 
C 

L"*"'"cOMMON  /  SSDATA  /  NSSPIS  ,  SDTIME  ,  SUBSTRC5000)  ,  SSTIME(5000) 
COMMON  /  REG  1 1  /  MAXPT  ,  MINPT  ,  VALMAX  ,  VALMIN 
LEVEL  2,  NSSPTS,  SOT  I  ME, SUbSTR, SST I  ME , MAXPT, Ml NPT, VALMAX, VALMIN 

L 

C  *  CALCULATE  MAXIMUM  AND  MINIMUM  VALUES 
C 

VALMAX  =  -99999 . 

VALMIN  =  99999. 

DO  9000  I  =  1» NSSPTS 

IF  (SUBSTK(I)  .GT.  VALMAX)  GO  TO  1000 
IF  (SUbSTR (I)  .LT.  VALMIN)  GO  1U  2000 
GO  TO  9000 
1000  CONTINOE 

*  SUBSTR ( 1 )  IS  NEW  MAXIMUM 

VALMAX  =  SUHSTR(I) 

MAXPT  =  1 
GO  TO  9000 
2000  CONTINUE 
C 

C  *  SUBSTR (I)  IS  NEW  MINIMUM 
C 

VALMIN  =  SUBSTR(I) 

MINPT  =  I 
9000  CONTINUE 
C 
C 

C  *  RETURN  MAXPT  =  INDEX  OF  MAXIMUM 

C  MINPT  =  INDEX  OF  MINIMUM 

C  VALMAX  =  MAXIMUM 

C  VALMIN  =  MINIMUM 

C 

9000  CONTINUE 
RETURN 
END 


PROGRAM  LISTING  (Continued) 
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APPENDIX  A. 2 

IMSL  ROUTINES  IN  IMSL2  LIBRARY 
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APPENDIX  A. 2 

IMSL  ROUTINES  IN  IMSL2  LIBRARY 


CALLED 

BY 

MOOSE 


CALLED  BY 
ANOTHER 
IMSL  ROUTINE 


RLFOR 

GFIT 

FTFREQ 

MDNOR 


RLFOTW 

RLDCW 

RLPOLY 

RLOPDC 

MDBETI 

MDBETA 

MLGAMA 

RLPRDI 

RLDCVA 

RLDOPM 

FTAUTO 

MDCH 

MGAMMA 

VERT ST 
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APPENDIX  B 

PROGRAM  INPUT 
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APPENDIX  B . 1 

JOB  CONTROL  LANGUAGE  FOR  RUNNING  THE  PROGRAM 
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APPENDIX  B.l 

JOB  CONTROL  LANGUAGE  FOR  RUNNING  THE  PROGRAM 

Job  Control  Language  for  running  MOOSE  UPDATE  file 
with  Level  2  IMSL  and  CALCOMP  Plotting  Routines. 

Job  name,Tl00 ,STMFZ . 

Account  card 
FILE , TAPE1 , RT=W . 

FILE ,OLDPL ,RT=S . 

REQUEST , TAPE10 , *PF . 

ATTACH, IMSL2L, ID  = 

LIBRARY, *,IMSL2L. 

BEGIN , ATTACH , PLOTLIB . 

COP YCR , INPUT , DEBG . 

REWIND, DEBG. 

GETPF ,OLDPL ,MOOSEOLDPL , ID  =  , ST=MFA. 

UPDATE, F,W,L=1. 

FTN , I , B=BMOOSE , L=0 , EL=F , D=DEBG , PL=6  5 00 . 

ATTACH ,TAPE1 , inductance  coil  data  file  name, ID  = 
MAP, OFF. 

B MOOSE . 

BEGIN , PLOT , CALCOMP , TAPE10 . 

EXIT. 

BEGIN , PLOT , CALCOMP , TAPE10 . 

EXIT. 

*EOR 

C#####DEBUG 

(Optional  user  debug  directives) 

*EOR 

*IDENT  TEMP 
*EOR 

user  input  data  deck 

*EOR 

*EOF 
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APPENDIX  B.2 

SUPPLEMENTAL  UPDATE  IDENT  FOR 
OBTAINING  PAGE  SIZE  PLOTS 
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APPENDIX  B . 2 

SUPPLEMENTAL  UPDATE  IDENT  FOR  OBTAINING  PAGE  SIZE  PLOTS 

*IDENT  PLT 
*D  PR-SCRIB4 .29,40 
XLNGTH=11 . 

YLNGTH=8 . 5 
XSIZE=8 .0 
YSIZEl=4. 

YSIZE2=1 . 5 
XBASE=2 . 

YBASE=3 . 0 

CALL  FIXSCA (XARRAY (1) ,NDPTS , XSIZE , XSCALE , XOR, XMAX, DX) 
CALL  FIXSCA (YARRYl(l) , NDPTS , YSIZEl , YSCALE , YOR, YMAX, DY) 
*D  PR-SCRIB4 . 50 

CALL  LABELA ( 0 . , DY , XOR , XMAX , YOR , YMAX , ONE , ONE ) 

*D  PR-SCRIB4 . 54 ,60 

CALL  FIXSCA (YARRY2 (1) ,NDPTS ,YS I ZE2 , YSCALE , YOR, YMAX, DY) 
YBASE=1 . 

*D  PR-SCRIB4 . 66 

CALL  PLT AXS ( DX , DY , XOR , XMAX , YOR , YMAX , MODE ) 

*D  PR-SCRIB4 . 86 

YMINl=YOR+ (YSCALE*6 . 3) 

*D  PR-SCRIB4 . 89 

CALL  PLTSYM  ( jZ(  •  15  ,  ISENT ,  0 .  ,XOR,YMINl) 

*D  PR-SCRIBL.  28 ,4j3 
XLNGTH=11 . 

YLNGTH=8 .5 
XSIZE=8 . 0 
YSIZE=6  ,)2f 
XBASE=2  .  jS 
YBASE=1 . 0 

CALL  FIXSCA ( XARRAY (1) , NDPTS , XSIZE , XSCALE , XOR, XMAX, DX) 
CALL  FIXSCA (YARRYl(l) , NDPTS , YSIZE , YSCALE , YOR, YMAX, DY) 
*D  PR-SCRIBL. 66, 69 
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SUPPLEMENTAL  UPDATE  IDENT  FOR  OBTAINING  PAGE  SIZE  PLOTS  (cont.) 

*D  PR-SCRIBL. 71 

YMAX1=YMAX+1 .  0*YSCALE 
*D  PR-SCRIBL. 74 

CALL  PLTSYM  (J? .  15 ,  ISENT , 0 .  ,XOR,YMAXl) 

*D  PR-SCRIBl .25,37 
XLNGTH=11. 

YLNGTH=8 . 5 
XSIZE=8 .0 
YSIZE=6.J3 
XBASE=2 . 0 
YBASE=1 . 0 

CALL  FI XS CA ( XARRAY ( 1 ) , NDPTS , XS I ZE , XSCALE , XOR , XMAX , DX ) 
CALL  FIXSCA (YARRAY (1) ,NDPTS ,YS I ZE , YSCALE , YOR, YMAX,DY) 

*D  PR-SCRIBl. 5 6, 5 9 
*D  PR-SCRIBl. 61 

YMAX= YMAX+ 1 . 0* YSCALE 
*D  PR-SCRIBl. 64 

CALL  PLTSYM  (J2I . 15  ,  ISENT,  )2f .  ,XOR,YMAXl) 

*D  PR-PAWS 2.256,268 
XLNGTH=11 . 

YLNGTH=8 . 5 
XSIZE=8 . 0 
YSIZE=6 . 0 
XBASE=2 . 0 
YBASE=1 . 0 

CALL  FIXSCA ( XARRAY ( 1 ) , NDPTS , XS I ZE , XSCALE , XOR , XMAX , DX) 
CALL  FIXSCA (YARRAY (1) , NDPTS ,YZISE, YSCALE, YOR, YMAX,DY) 

*D  PR-PAWS 2.287,29)21 
*D  PR-PAWS2 . 292 

YMAX1=YMAX+1  .)2(*YSCALE 
*D  PR-PAWS2 .295 

CALL  PLTSYM()2f.  15,  ISENT, 0.  , XOR, YMAX1) 
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APPENDIX  B. 3 

INPUT  DATA  DESCRIPTION 
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APPENDIX  B  .‘3 
INPUT  DATA  DESCRIPTION 


Input 

Level  Subroutine/Overlay  Card  Format  Variable 


Description 


IAINPT  (A, 0,0) 


USRNAM  (2) 


Called  only  once* 
User  name,  2  names 
each  with  a  maximum 
of  10  characters. 


2 

RDPURP  (A, 0,0) 

1 

15 

2 

NLINES+1 

6A10 

Called  only  once. 

NLINES  Number  of  cards  used  to 

enter  purpose  of  run. 
PURPOS(6)  Purpose  of  run  -line  1 
Purpose  of  run  -line 
NLINE 


3  ANALIZ  (A, 0,0)  Repetitively  used  to 

control  analysis. 

1  12  IANOPT  Analysis  option 

=0  No  further  analysis. 

Program  terminates. 
=1  User  controlled 
analysis . 

Control  passes  to 
subroutine  KITTY8. 
=2  Complete  data  record 
analysis . 

Control  passes  to 
subroutine  COMPLT. 


****  END  OF  INPUT  DECK  IF  IANOPT  =  0  **** 
****  NEXT  CARD  IS  INPUT  LEVEL  4  OF  **** 
****  KITTY 8  OR  COMPLT  INPUT  OPTIONS  **** 
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INPUT  DATA  DESCRIPTION  (cont.) 


Input 

Level 

4 


Subroutine/Overlay  Card  Format  Variable  _ Description 

KITTY8  (A, 0,0)  Optionally  called  by 

ANALIZ 
IANOPT  =  1 

1  12  SUBOPT  Analysis  option 

=0  No  further  analysis. 
Control  returns  to 
ANALIZ. 

=1  Select  data  subset. 
Control  passes  to 
subroutine  SUBSET. 

=2  Power  spectrum. 

Control  passes  to 
subroutine  PAWS. 

SUBOP T=1  must  be 
called  before  this 
option. 

=3  Forward-backward  fits. 
Control  passes  to 
subroutine  GIRAFE. 

Not  currently  avail¬ 
able. 

=4  Polynomial  and  per¬ 
iodic  fits,  noise 
analysis . 

Control  passes  to 
subroutine  THEFFT. 
SUBOPT=l  must  be 
called  before  this 
option. 

SUBOPT=ll  must  be 
called  before  this 
option. 

=5  Smoothing. 

Control  passes  to 
subroutine  SMOOTH. 
SUBOPT=l  must  be 
called  before  this 
option. 

=11  First  and  second 
derivatives  of 
entire  data  set. 
Control  passes  to 
subroutine  DERIV. 

=12  Detection  of  two 
epochs. 

Control  passes  to 
subroutine  DETECT. 
SUBOPT=ll  must  be 
called  before  this 
option. 
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INPUT  DATA  DESCRIPTION  (cont.) 

Input 

Level  S  ub  rout ine/O ve  r 1 ay  Card  Format  Variable  Description 

****  END  OF  KITTY8  INPUT  IF  SUBOPT  =  0  **** 

****  NEXT  CARD  IS  INPUT  LEVEL  5  **** 

****  SEE  APPROPRIATE  SUBROUTINE  BELOW  **** 


5  SUBSET  (C,2,0) 
SUBPAR  (C, 2,0) 


5  PAWS  (0,3,0) 


5  GIRAFE  (E , 4 , 0) 


1 

F10.4 

TINIT 

F10.4 

TLAST 

110 

NODTS 

1  12  PLOTOP 


12  LAGOPT 


Optionally  called  by 
KITTY 8  with  SUBOPT=l 
Initial  time  of  data 
subset,  milliseconds . 
Final  time  of  data 
subset,  milliseconds . 
Number  of  2.5  micro¬ 
second  intervals  be¬ 
tween  points  in  data 
subset. 


Optionally  called  by 
KITTY8  with  SUBOPT=2 

Plotting  option  for 
power  spectrum. 

=0  Plot  off. 

=1  Plot  on. 

Frequency  option. 

=0  Frequencies  resolu¬ 
tion  is  maximum 
permitted,  multiples 
of  .25,  .5,  1.0  or 
2.0  kHz. 

=1  Frequencies  are  mul¬ 
tiples  of  1  kHz. 


Optionally  called  by 
KITTY8  with  SUB0PT=3 
Not  currently  avail¬ 
able. 

No  input  required. 


5  THEFFT  (F,5,0) 


Optionally  called  by 
KITTY8  with  SUBOPT=4 

I  12  PLOTOP (1)  Plotting  option  for  data 

fits. 

=0  Plot  off. 

=1  Plot  on. 
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INPUT  DATA  DESCRIPTION  (cont.) 


Input 

Level  Subroutine/Overlay  Card  Format  Variable  _ Description 


PAWS  (1,5,2) 


12 

PLOTOP (2) 

Plotting  option  for 
entire  fitted  data  set 
=0  Plot  off. 

=1  Plot  on. 

Called  by  THEFFT 

2 

12 

PLOTOP 

Plotting  option  for 
power  spectrum. 

=0  Plot  off. 

=1  Plot  on. 

5  SMOOTH  (G,6,0) 


1  12  SMOPT 


Optionally  called  by 
KITTY 8  with  SUB0PT=5 
Smoothing  option. 

=1  11-point  local 
median . 

=2  11-point  local  mean. 


5  DERIV  (L, 7 ,0) 

1  12 

5  DETECT  (N, 10 , 0) 


PLOTOP 

Optionally  called  by 
KITTY8  with  SUBOP T= 11 
Plotting  option  for 
first  derivatives. 

=0  Plot  off. 

=1  Plot  on. 

Optionally  called  by 
KITTY8  with  SUBOPT=12 
No  input  required. 

****  NEXT  CARD  IS  INPUT  LEVEL  4  FOR  NEXT  **** 
****  ANALYSIS  OPTION  OF  SUBROUTINE  KITTY8  **** 


4  COMPLT  (A, 0,0) 

DERIV  (L, 7,0)  1 


Optionally  called  by 
ANALIZ 
IANOPT  =  2 

12  PLOTOP  Plotting  option  for 

derivatives. 

=0  Plot  off. 

=1  Plot  on. 
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INPUT  DATA  DESCRIPTION  (cont.) 


Input 

Level 

Subroutine/Overlay 

Card 

Format 

Variable 

Description 

COMPLT  (A, 0,0) 

2 

12 

IDIV 

Option  for  division  of 

Regime  I 

=“1  Divide  Regime  I 
=0  Do  not  divide 


Regime  I 

=1  No  further  analysis 
Control  returns  to 
ANAL I Z 


**** 

NEXT 

CARD 

IS 

INPUT  LEVEL  3  IF  IDIV=1 

**** 

**** 

NEXT 

CARD 

IS 

CARD  4  IF 

IDIV=0 

**** 

**** 

NEXT 

CARD 

IS 

CARD  3  IF 

IDIV=-1 

**** 

COMPLT  (A, 0,0) 

3 

F8.4 

SUBTIM 

THEFFT  (F, 5 , 0) 

4 

12 

PLOTOP  (1) 

12 

PLOTOP  (2 ) 

PAWS  (1,5,2) 

5 

12 

PLOTOP 

PAWS  (D, 3,0) 

6 

12 

PLOTOP 

12 

LAGOPT 

****  REPEAT  CARDS  4,  5,  AND  6  FOR  **** 


Time  of  subdivision  of 
Regime  I  in  milliseconds 

Plotting  option  for 
Regime  I  data  fits. 

=0  Plot  off. 

=1  Plot  on. 

=0  Required  input. 

Plotting  option  for 
Regime  I  power  spectrum. 

=0  Plot  off. 

=1  Plot  on. 

Plotting  option  for 
Regime  I  noise  power 
spectrum. 

=0  Plot  off. 

=1  Plot  on. 

Frequency  option  for 
Regime  I  noise. 

=0  Frequencies  resolu- 
tion  is  minimum 
permitted,  multiples 
of  .25,  .5,  1.0  or 
2 . 0  kHz 

=1  Frequencies  are 

multiples  of  1.0  kHz. 


****  PART  2  OF  REGIME  I  IF  IDIV=-1  **** 
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INPUT  DATA  DESCRIPTION  (cont.) 


Input 

Level 

Subroutine/Overlay 

Card 

Format 

Variable 

Description 

PAWS  (D, 3,0) 

7 

12 

PLOTOP 

Plotting  option  for 

Regime  II  power 
spectrum. 

=0  Plot  off. 

=1  Plot  on. 

12 

LAGOPT 

=1  Required  input. 

THEFFT  (F, 5 ,0) 

8 

12 

PLOTOP (1) 

Plotting  option  for 
Regime  III  data  fits. 
=0  Plot  off. 

=1  Plot  on. 

12 

PLOTOP (2) 

Plotting  option  for 
entire  fitted  data  set 
=0  Plot  off. 

=1  Plot  on. 

PAWS  (1,5,2) 

9 

12 

PLOTOP 

Plotting  option  for 

Regime  III  power 
spectrum. 

=0  Plot  off. 

=1  Plot  on. 

PAWS  (D, 3,0) 

10 

12 

PLOTOP 

Plotting  option  for 

Regime  III  noise 
power  spectrum. 

=0  Plot  off. 

=1  Plot  on. 

12  LAGOPT  Frequency  option  for 

Regime  III  noise. 

=0  Frequencies  reso¬ 
lution  is  maximum 
permitted,  multiples 
of  .25,  .5,  1.0  or 
2.0  kHz. 

=1  Frequencies  are  mul¬ 
tiples  of  1.0  kHz. 

****  NEXT  CARD  IS  INPUT  LEVEL  3  **** 


*  *  *  * 


ANALYSIS  OPTION  OF  SUBROUTINE  ANALIZ  **** 


APPENDIX  B . 4 

SAMPLE  INPUT  DATA 


119 


APPENDIX  B.4. 
SAMPLE  INPUT  DATA 


This  data  deck  will  produce  the  output  found  in 
Appendix  C.2,  Sample  Output. 


*EOR 

"SAMPLE" 

1 

TEST  -  CHECK  COMPLT  LEG 
2 
1 
0 

1  0 
1 

1  0 
1  1 
1  1 
1 

1  0 

*EOR 
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APPENDIX  C 

PROGRAM  OUTPUT 
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APPENDIX  C.l 

OUTPUT  DATA  DESCRIPTION 
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APPENDIX  C.l 
OUTPUT  DATA  DESCRIPTION 


User ' s  Name 
Purpose  of  Run 
Analysis  Option 
Data  Record  Description 
start  time 
number  of  point 
point  interval 
Data  Record  Classification 
saturated  points 
classification 

Portions  of  a  Complete  Analysis  Follows,  KITTY8  User's 
Option  will  vary 

Derivative 

no  printout 

Epoch  Detection 

index,  time,  values  of  Epochs  1  and  2 
Optional  Division  of  Regime  I 

Analysis  of  Regime  I,  III,  or  Parts  of  Regime  I 

Regime  Fitting  and  Analysis  (THEFFT  Subroutine) 

polynomial  coefficients  for  x,  x2,  x3,  x4,  x° 
Power  Spectrum  (PAWS  Routine) 

mean  and  variance  of  waveform 
lag  and  autocovariance 
frequency  and  power 
ten  predominant  frequencies 
Periodic  Curve 

frequency,  cosine  and  sine  coefficients 
Fitted  Curves 

index,  time,  data,  polynomial  fit,  data  minus 
polynomial  fit,  periodic  fit,  noise 

Analysis  of  Noise 

mean,  standard  deviation 

X2  test,  x2  ami  Q  statistics 
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OUTPUT  DATA  DESCRIPTION  (cont.) 


Shock  Detection 

for  points  at  least  2.5  standard  deviations  from 
the  mean:  index,  time,  noise  value,  standard 
deviations  from  mean,  a  if  more  than  3., 

standard  deviations  from  mean  of  derivative 

Analysis  of  Regime  II 

maximum  and  minimum 
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APPENDIX  C.2 

SAMPLE  OUTPUT 
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APPENDIX  C. 2 
SAMPLE  OUTPUT 


The  following  plots  were  produced  using  the  input 
data  listed  in  Appendix  B.4,  and  the  UPDATE  deck  for  obtaining 
page  size  plots,  in  Appendix  B.2. 

1.  Top:  Plot  of  data  as  stored  on  tape. 

Bottom:  Plot  of  first  derivative  of  data. 

2.  Power  spectrum  of  Regime  I. 

3.  Solid  line:  Data  on  Regime  I. 

th 

Broken  line:  41'  degree  polynomial  fit  to  data  in 
Regime  I . 

4.  Solid  line:  Regime  I  data  after  subtraction  of 

polynomial . 

Broken  line:  Fourier  series  fit  using  ten  predominant 
frequencies . 

5.  Regime  I  noise;  data  with  polynomial  and  Fourier  fit 
subtracted. 

6.  Power  spectrum  of  Regime  I  data  after  subtraction  of 
polynomial . 

7.  Power  spectrum  of  Regime  II. 

8.  Power  spectrum  of  Regime  III. 

9.  Solid  line:  Data  in  Regime  III. 

j 

Broken  line:  4r  degree  polynomial  fit  to  data  in 
Regime  III. 

10.  Solid  line:  Regime  III  data  after  subtraction  of 

polynomial . 

Broken  line:  Fourier  series  fit  using  ten  predominant 
frequencies. 

11.  The  Regime  III  noise;  data  after  subtraction  of 
polynomial  and  Fourier  fit. 

12.  Entire  data  set  with  noise  from  Regimes  I  and  III 
removed . 

13.  Power  spectrum  of  Regime  III  after  subtraction  of 
polynomial . 
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50.5  50.6  50.7  50-8  50. 


50.4  50.5  50.6  50.7  50.8  50. 
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53.4  53.8  54. 


53.0  53.4  53.8  54.2  54.6  55-0  55 


APPENDIX  C . 3 


EXECUTION  ERROR  CODE  DESCRIPTION 


APPENDIX  C . 3 

EXECUTION  ERROR  CODE  DESCRIPTION 

PEKROR 

ERROR 

CODE  SUBROUTINE/OVERLAY  ERROR  DESCRIPTION 


5  ANALIZ  (A, 0,0)  IMPROPER  ANALYSIS  OPTION. 

IANOPT  <  0 


8  ANALIZ  (A, 0,0)  IMPROPER  ANALYSIS  OPTION. 

IANOPT  >  2 


51  SUBPAR  (C , 2 , 0)  IMPROPER  INPUT  TO  SUBROUTINE. 

TLAST  <  TINIT,  DATA  SUBSET  AND  POINTS 


52  SUBPAR  (C,2 ,0)  IMPROPER  INPUT  TO  SUBROUTINE. 

POINT  INTERVAL  BETWEEN  DATA  SUBSET  POINTS, 
NODTS  <  0 


53  SUBSET  (C,2 ,0)  CALCULATION  ERROR. 

IS  <  0,  DATA  SUBSET  STARTING  INDEX 


54  SUBSET  (C,2 ,0)  CALCULATION  ERROR. 

IE  <  IS,  ENDPOINT  INDICES  OF  DATA  SUBSET 


55  SUBSET  (C,2 ,0)  DATA  SUBSET  ERROR. 

NSSPTS  <  0,  NUMBER  OF  SUBSET  DATA  POINTS 


57  PAWS  (D , 3 , 0)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ. 

IND (2)  <  3,  NSSPTS  <  3 


58  PAWS  (D, 3 ,0)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ. 

IND (4)  <  2,  M  <  2 


59  PAWS  (D , 3 ,0)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ. 

IND (4)  >  IND (2),  M  >  NSSPTS 
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EXECUTION  ERROR  CODE  DESCRIPTION  (cont.) 


PERROR 

ERROR 

CODE  SUBROUTINE/OVERLAY  ERROR  DESCRIPTION 


60  SUBSET  (C,2,0)  CALCULATION  ERROR  IN  DATA  SUBSET  INITIAL 

POINT.  CALCULATED  POINT  DOES  NOT  CORRESPOND 
TO  ACTUAL  INITIAL  POINT 


61  SUBSET  (C,2 ,0)  CALCULATION  ERROR  IN  DATA  SUBSET  FINAL 

POINT.  CALCULATED  POINT  DOES  NOT  CORRESPOND 
TO  ACTUAL  FINAL  POINT. 


63  THEFFT  (F,5,0)  UNABLE  TO  REPLACE  DATA  INTO  DATA  RECORD. 

SDTIME  ?  DTIME 


64  THEFFT  (F,5,0)  CALCULATION  ERROR  IN  DATA  SUBSET  INITIAL 

POINT.  CALCULATED  POINT  DOES  NOT  CORRESPOND 
WITH  ACTUAL  INITIAL  POINT. 


71  COMPLT  (A, 0,0)  IMPROPER  REGIME  I  DIVISION  OPTION. 

IDIV  J-  -1,  0,  1 


73  COMPLT  (A, 0,0)  IMPROPER  DIVISION  TIME  OF  REGIME  I. 

SUBTIM  NOT  IN  REGIME  I. 


79  KITTY8  (A, 0,0)  IMPROPER  OPTION  FOR  ANALYSIS  OF  DATA 

RECORD  OR  DATA  SUBSET.  SUBOPT  <  0 
SUBOPT  1  1,  2,  3,  4,  5,  11,  OR  12 


81  THEFFT  (F,5,0)  ERROR  IN  ARRAY  DIMENSION  DATA  SUBSET. 

ARRAY  DIMENSION  TOO  LARGE 


82  SCRIB1  (K, 3, 1)  ERROR  IN  ARRAY  DIMENSION. 

NDPTS  >  100,  TOO  MANY  POINTS  TO  PLOT 


83  SCRIBL  ( J,5 , 3)  ERROR  IN  INPUT  TO  PLOTTING  ROUTINE. 

NARRAY  <  1  OR  NARRAY  >  3,  NUMBER  OF 
DATA  ARRAYS  TO  PLOT  ON  A  GRAPH 
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EXECUTION  ERROR  CODE  DESCRIPTION  (cont.) 


PERROR 

ERROR 

CODE  SUBROUTINE/OVERLAY  ERROR  DESCRIPTION 


84  SCRIBL  ( J, 5 , 3)  ERROR  IN  ARRAY  DIMENSION.  NDPTS  >5000, 

TOO  MANY  POINTS  TO  PLOT 


85  THEFFT  (F,5,0)  EXECUTION  ERROR  IN  IMSL  SUBROUTINE  GIFT. 

IER  >  0,  SEE  IMSL  MANUAL. 


87  POLYOR  (H,5 ,1)  EXECUTION  ERROR  IN  IMSL  SUBROUTINE  RLFOR 

IER  =  33,  SEE  IMSL  MANUAL 


88  POLYOR  (H, 5,1)  EXECUTION  ERROR  IN  IMSL  SUBROUTINE  RLFOR 

IER  >  0,  SEE  IMSL  MANUAL 


90  PAWS  (1,5,2)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ . 

IND (2)  <  3,  NSSPTS  <  3 


91  PAWS  (1,5,2)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ. 

IND (4)  <2,  M  <  2 


92  PAWS  (1,5,2)  ERROR  IN  INPUT  PARAMETER  ARRAY  TO  IMSL 

SUBROUTINE  FTFREQ. 

IND (4)  >  IND (2) ,  M  >  NSSPTS 


93  SCRIB2  (1,5,2)  ERROR  IN  ARRAY  DIMENSION. 

NDPTS  >  100,  TOO  MANY  POINTS  TO  PLOT 


94  SMOOTH  (G, 6,0)  IMPROPER  SMOOTHING  OPTION. 

SMOPT  <  1,  OR  SMOPT  >  2 


95  SCRIB4  (M, 7,1)  ERROR  IN  ARRAY  DIMENSION. 

NDPTS  >  5000,  TOO  MANY  POINTS  TO  PLOT 
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